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I. INTRODUCTION 



A. EACKGECONE 



The Fcurier transform has been with us many years and 
the unigueness of the transform allows Fourier analysis 
techniques developed in one field to be applied to many 
diverse areas. Typical application areas of the Fcurier 
transform include: 



Linear Systems - The Fourier transform of the 
output cf a linear system is given by the product 
cf the system transfer function and the Fourier 
trarsfcrm of the input signal. 



Antennas — The field pattern cf an antenna is 
given by the Fourier transform of the antenna 
current illumination. 



Optics - Optical systems have the property that a 
Fourier transform relation exists between the 
light amplitude distribution at the front and back 
focal planes cf a converging lens. 



Random Process - The power density spectrum of a 
random process is given by the Fcurier transform 
cf the auto-ccrrelation function cf the process. 



Probability - The characteristic function of a 
randcm variable is defined as the Fourier 
transform of the probability density function of 
the randcm variable. 



Quantum Physics — The "uncertainty principle" in 
quantum theory is fundamentally associated with 
the Fcurier transform since particle momentum and 
pcsiticr are essentially related through the 
rcurier transform- 



Eouncary- Value Problems - The solution of partial 
differential equations can he obtained by means of 
the Fcurier transform. 
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Although Fourier analysis allows one to more easily 
examine a function from another point of view,i.e., the 
transfom domain, the methods available before 1965 to 
confute the Fcurier coefficients were very time consuming 
and costly. Then in 1965 Cooley and Tukey published their 
mathematical algorithm which computed the Fourier 
coexficents with much less computation effort than had teen 
required in the past. This method has become known as the 
"fast Fourier transform 11 or FFT and has produced major 
changes in computational techniques used in digital spectral 
analysis, filter simulation and related fields. 

Although the FFT is the most well known algorithm to 
compute the Fourier coefficients, it is net the only method 
for computing the "discrete Fourier transform" (DFT) . The 
chirp-z-transf erm (CZT) is an algorithm for evaluating the 
z-transferm cf a finite duration sequence along certain 
general contours in the z-plane. Evaluating the DFT cf the 
seguence is a special case of the CZT, i.e., the z-transform 
cf the finite duration sequence is evaluated on the unit 
circle in the z-plane. Although the CZT is not quite as 
efficient as the FFT, it eliminates many of the restrictions 
of the latter. These restrictions are listed at the end of 
the CZT algorithm derivation. The "prime transform" is an 
algorithm that fellows directly from the CZT derivation. It 
is based upon "N", the number cf samples cf data, being an 
odd prime. 

The CZT and prime transform algorithms are receiving 
considerable attention for applications in spectral 
analysis . This has been brought about because the bulk of 
the computations of these algorithms is performed by a 
transversal filter, which is easily implemented by charge 
transfer devices (CID) or surface acoustic wave (SAW) 
devices. Hardware and software implementations of the FFT, 
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CZT, and prime transform algorithms ar€ discussed later in 
this chapter. 

B. BASIC CCKCEETS OF FOURIER TRANSFORM ANALYSIS 

The esserce of the Fourier transform of a waveform is to 
decompose or separate the waveform into a sum of sinusoids 
of different frequencies. The pictorial representation of 
the Fourier transform is a diagram which displays the 
amplitude anc frequency of each of the determined sinusoids. 

Mathematically, the AMELITUDE SPECTRAL CENSITY or more 
generally the Fourier transform of x(t) is given by the 
relationship 



where x (t) is the waveform to be decomposed into a sum of 



The INVERSE FOURIER TRANSFORM is given by the relationship 



The Fourier transform and its inverse may conveniently be 
written as 



oo 




( 1 - 1 ) 




oo 




v(t) <3=j> v(f) 

FOURIER TRANSFORM E AIR 
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For discrete signal processing the equations take cn the 
following fcim. For a periodic function x(t) the Fourier 
transfer! is defined by 



X(k) 



N-1 




(2tt/N) nk 



n = 0 



, k=0 , 1 . ,N- 1 (1-3) 



The "true" freguency component is given by 



X (k/NT) 



N-1 

•£ 

n = 0 



(nT) e 



- j (2 tt/N) nk 



(1-4) 



Egs . [ 1-3 ] and [1-4 ] are known as the EISCEETE FCU5IEB 
TE2NSFQES (EM). The inverse EFT is given by the following 
relation 



N-1 



Y + j (2tt/N) nk 
X(k)e 



k=0 



,n=0, 1,...,N-1 (1-5) 



ard for "true" time 



N-1 

x(nl) = (1/ip) X (k/Hl) e 

k=0 



+ j (2 tt/N) nk 



( 1 - 6 ) 
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where 



N — > total number of samples cf time 
and samples of frequency 

T — > sample period 

n — > sample period index 

k — > harmonic index 

P = ( 1/NT) fundamental frequency 
c 

Fcr future reference, the basic properties cf the 
continuous Fourier transform and discrete Fourier transform 
are sumnarimed in Table-I. Fig.[1] gives a pictorial 
representation of some of the most commonly used Fourier 
transform pairs. 
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Frequency domain 
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Figure 1 - FOURIER TRANSFORM PAIRS 
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Figure 1 - CONTINUED 



Time domain Frequency domain 
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Figure 1 - CONTINUED 



C. IKPLEfEMATIGN PRCELEHS 



lh€ twc problems most often encountered in using the 
discrete fctrier transform are aliasing and leakage. 



1 . ^ lias i n_g 



Aliasing refers to the fact that the high-f reg uency 
components cf a time function can impersonate low 
freguencies if the sampling rate is too low. This problem 
is corrected by sampling the signal at a rate > the Nyguist 
Freguency, i.e., demanding that the sampling rate be high 
enough fcr the highest freguency present to be sampled at 
least twice during the period. Fig.[2] gives a pictorial 
representaticn of a signal without aliasing. 




(e) 



f\ 



T7 



4h(t) &(t) 



-V 



r 




figure 2 - Fourier Transform of a Waveform Sampled at the 
Nyguist Sampling Rate. 
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ih€D the sampling period T of Fig.[2-c] is increased 
as shown in Fig.[3-c0, the equidistant impulses of (f) 
become acre closely spaced. Fig. [3-d]. Eecause of the 
decreased spacing of the freguency iapulses, their 
convolution with the frequency function, a (f ) , Fig.[3-b 
results in the overlapping waveform illustrated results in 
the overlapping waveform illustrated in Fig.[3-f]. This 
overlap is the result of aliasing by the high frequency 
ccapcnen ts . 



The aliasing phenomenon will occur anytime a sampled 

analog signal contains freguency components greater than 

half the saapling freguency, (f ) /2. The aliasing frequency 

s 

relationship is such that all frequencies higher than the 

analysis bandwidth [f = (f )/2] "fold" down into that 

max s 

band in an accordian-like manner with a fold at every 



multiple cf f 




(e) 



h(t)£<t) 



1 A 




V "T 

Eigure 3 - Fourier Transform of a Wavefcra Sampled at Less 
Than the Nyguist Sampling Fate. 



leakace 



If a periodic , band limited function is sampled and 
truncated tc consist of other than an integer multiple of 
the period, the resulting discrete and continuous Fourier 
transforms %ill differ considerably. This problem is 
inherent in the Fourier analysis of any finite record of 
data. The record has been formed by looking at the actual 
signal for T seconds and by neglecting everything that 
happened before or after this period. As shown in 
Fig.[4-cj, this is eguivalent to multiplying the signal by a 
rectangular window. 

If the continuous Fourier transform of the pure 
cosine wave of Fig.[4-a] had been found, its contribution 
wculd have teen limited to a pair of impulse functions on 
the frequency axis, 5ig.[4-a]. 
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figure 4 - DFT of a Eand-limited Periodic 
Waveform :Truncation Equal tc an Integer Multiple. 
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The effect cf truncation at other than a multiple of 
the period cf the sampled signal is to create a periodic 
function with sharp discontinuities as shewn in Fig. [5]. 
The multiplication by the data window in the time domain is 
eguivalent tc performing convolution in the freguency 
domain. Consequently, the frequency function is no longer a 
single impulse but rather a continuous function of freguency 
with a local maximum centered at the original impulse and a 
series of sptiious peaks termed sidelobes. These sidelcbes 
are respcDsitle for the additional frequency components 
"leakage", which cccur after frequency domain sampling. An 
expanded view cf sampling in the freguency domain is shown 
in Fig .[ 6 ]. 
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figure 5 - EFT cf a Periodic Waveform: Truncation Interval 
Not Equal to an Integer Multiple 
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figure 6 - Expanded Illustration cf the Convolution of 
Figure [ 5e ] . 



The csual approach tc reduce the ancunt of “leakage” 
through these sidelohes consists cf applying a data window 
or weighting function to the time series, which has lower 
sidelohes in the frequency domain than the rectangular data 
window. The process of weighting consists of 
amp lit ud e-mo culati rg the signal prior to Fourier analysis. 
The weighting process is illustrated in Fig. [7], 



TIME 

FUNCTION ^ 


FOURIER 

TRANSFORM 


FOURIER TRANSFORM 
OF PRODUCT EQUALS 


w <y * 


TRANSFORM OF THE 
TIME FUNCTION 



CONVOLVED WITH 



WEIGHTING 

FUNCTION 



TRANSFORM OF THE 
WEIGHTING FUNCTION 



figure 7 - The Weighting Function Process. 
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lhere are many such windows but one of the simplist 
and most often used is the Hanning function. This function 
is a cosine-bell curve given by 

w (t) = 1 - 1ccs 2^t 0<t<Tc ( 1-7) 

where T^ is the truncation interval. The magnitude of the 
Fourier transform cf the Hanning function is given by 



|W(I) I = 1Q(f) ♦ 1[Q(f+1 ) + C(f-1 ) ] (1-8) 

^ 3 Tc Tc 

where C(f) = [ sin (rrT c f ) ]/TTf 

The Hannirg function Fourier transform fair is shewn in 
Fic.(e]. 

Windowing of the input frame by the Hanning function 
reduces "leakage" by imposing a guasi-pericdicit y or the 
input sicral and eliminating the discontinuities between 
periods. The Hanning function widens the main lobe of the 
(sin x) /x shape and reduces the spectral amplitude by a 
little more than 4dE, but it causes the "interfering" 
sidelcbes to fall off at 18d3 per octave instead cf 6dE as 
in the rectangular window. 
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AMPLITUDE 





Figure 



6 



Eanning Function Fourier Transform Pair 



It should be remembered that the non-zero frequency 
components are considerably EBOADENED cr SMEABED with 
respect tc the impulse function. In general, the mere the 
leakage is reduced, the broader or more smeared the results 
of the discrete Ecurier transform appear. Iable-II contains 
seme of the acre ccmmcnly used data windows. 
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TABLE-II 
Data Windows 



WEIGHTING FUNCTION 


3 -dB 

BW 


NOISE 

BW 


HIGHEST 
SIDELOBE 
LEVEL (dB) 


ASYMPTOTIC 

ROLL-OFF 

(dB/OCTAVE) 


RECTANG 


JLAR 


1/B 


0 . 85 B 


IB 


-13 


6 


TRIANGL 


E 


1.25B 


1.35B 


-26 


12 


COSINE 




1.17B 


1 . 26 B 


-23 


12 


HANNING 


= (COSINE)^ 


l.UB 


1.5B 


-32 


18 


(COSIIiE 


)3 


1.613 


1.73B 


-39 


24 


(COSUE 




1.88B 


1.9B 


-48 


30 


HAVING 


= (COSINE ) 2 + 

3/6 PEDESTAL 

\ — ? 


1.3B 


1 . 36 B 


-42 


6 dB/OCT 

BEYOND 

5B 
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SPEC! BUM ANALYZES PERFORMANCE FIGURES OF MERIT 



1 • fa spl inq Frequ ency and Anti-Aliasing Fi lte r 

The sampling frequency must satisfy the Nyguist 
sampling theorem or aliasing will result. For data that has 
conpcnents beyond the analysis range, the components are 
attenuated ty means of a lowpass filter, known as an 
"anti-aliasing" filter. A typical "fall off" rate for 
anti-aliasing filters in today*s equipment is 120 dE/octave. 
In practice, a sampling rate of approximately three times 
the width cf the desired analysis range is a good 
compromise. Sampling at higher rates requires more 
equipment, which is undersirafcle fer eccncoiic reasons. 
Sampling at lewer rates fails, to take into account the 
actual fall-cff of available anti-aliasing filters, anc thus 
dees not eliminate aliasing completely. 

2 . Bandwidth (Freq uen cy Rangef 

The bandwidth determining parameters of a spectrum 
analyzer are: 

1. The memory capacity (which determines the 
cumber of samples available fer processing) . 

2. The sampling frequency (which, in combination 

vitfc determines the length cf signal (ir 

seconds) which is available fer processing). 

3. The shape of the "window function" also known 
as the "weighting function". 
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3 



. Egscluticn ce 3-dB Bandwidth 

In Table-II the 3-dB bandwidth for various weighting 
functions were compared using the relation 

B (H 2) = 1 

Reighting-f unct icn length (sec) 

The 3-dE bandwidth is a good measure cf a spectrum 
analyzer's ability to resolve two equal-amplitude adjacent 
sire waves. Studying Table-II, it appears that the 
rectangular weighting function has the best resolution but 
in reality, the response characteristic is limited in 
usefulness by the side-lobe structure. In practice, it is 
mere important to resolve unegual than egual amplitude 
adjacent frequency components. 

4 • ]jcis e Bandwidth 



Sc iie-times called "effective bandwidth", it is used 
to convert cr normalize power spectrum measurements to power 
spectral density (pewer per unit frequency) . The ncise 
bandwidth is the bandwidtn cf a hypothetical rectangular 
filter, which passes a signal with the same mean-sguare 
value as the actual filter when the filter input is white 
Gaussian rcise. 

5 • highest Si delobe Level and A sv apt o tic R cll -of f 

It is useful to specify these parameter, since it 
provides seme indication of the prominence cf the side-lobe 
structure. The largest sidelobe is generally the first cne. 
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with the Hamming function being among the exceptions. The 
third side lcte of the Hamming function is the largest. 

The presence of side lobes makes it difficult to 
distinguish between a high-level frequency component and a 
lew-level component which is close in frequency. For this 
reason, rectangular weighting, which leads to a spectrum 
response egual to the true Fourier transform of the 
tin e-tru nca ted signal, is used only occasionally when 
analyzing trarsients or shocks. 

6 . f icgle^Sid e d S pect ra 

Although it is net an actual measure of performace, 
a discussion of single-sided spectra is warranted at this 
time so that the results obtained frem the various 
algorithms is fully understood. 

The firite discrete transform presumes that the 
input waveform is the sum of cosine-phase and sine-phase 
freguency components, all harmcnically-related to the lowest 
freguency that can be recognized. That fundamental is the 
freguency that completes one cycle in the length of the 
input frame. The calculated harmonics range from <d-c) 
through 1 (the fundamental) , 2,3, etc. , to N/2, the maximum 
recognizable tarmcnic, which completes one cycle in the two 
intervals associated with two pents of the input (i.e., two 
sa ii pies per cycle) . 

Ir the general case, the discrete Fcurier transform 
can operate on a complex time series input function and 
produce a ccnplex spectrum. In such a case, an input 
furcticn cf N complex values (i.e., 2N independent real and 
imaginary parts) yields a spectrum of 2N independent real 
and imaginary parts representing N complex frequency 



35 



coefficients centered about 0 (d-c) and ranging from harmonic 
-N/2 to +N/2. The "negative- frequency H components are an 
integral and necessary part cf the function. 

A real input function, as an electrical waveform, 
may te considered for this general case as a complex 
function whose imaginary parts are all zero. The discrete 
transform cf such an input yields a complex spectrum as 
abcve. For this special-case input however, the spectrum 
will be symmetrical atcut the center or d-c point. The 
nega tive-f recuency parts from 0 to -N/2 are an image cf the 
positive f r e guenci es-, the real or cosine parts having even 
symmetry arc the imaginary, sine parts odd. The 

necative-fr ecuency half thus becomes redundant, and an input 
of N real points yields N independent real and imaginary 
parts cf S/2 complex spectral coefficients. 

Scnet heless , the negative-frequency coefficients are 

still a necessary part of the spectrum in that they 

represent half the amplitude of the harmonic components. In 

ether words, if the input function contains a cosine-phase 

ccmpcnent cf amplitude A at frequency 3f , then the 

c 0 

resultant pc sitive-f requency coefficient for frequency 3f^ 
will have amplitude A/2 and the n egative- f recuency 

th 

coefficient fer -3f will also have amplitude A/2. The 0 
0 

or d-c coefficient, however, dividing or being "shared" 

between the positive and negative halves cf the spectrum, 
will have full amplitude A for an input d-c level of A, 
which is the point cf this discussion. 
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E. 5 FEC1EAL A NALY SIS ALGORITHMS 



1 • 11 e Fast Feu ri er Transfora 



The fast Fourier transform (FFT) is a highly 



efficient prccedure for computing the DFT cf a time series. 
It takes advantage cf the fact that the calculation of the 
coefficients cf the DFT can be carried cut iteratively, 
which results in a considerable savings of ccmputaticn time. 
Tc gain insight to the FFT algorithm, a general developement 
is presented fcllcwed by several examples. A more in depth 
study of the FFT may be found in Refs. [6 ], [16], and [17], 



a 



General Developement cf the FFT Algorithm 



The discrete Fourier transform Eg. [1-3] is 



N - 1 





n = 0 



which maj he written as 



N-1 




(1-9) 



n =0 
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where 



M 



-j2tt/N 

€ 



( 1 * 10 ) 



i . 6 



nk 

Is can be seen, W is periodic and of period N; 



(n+mti) (k+lN) nk 

» = W , m,l= 0,±1 , (1*11) 



In carrying cut the FFT algorithm, the symmetry and 
nk 

periodicity cf Vi are exploited to achieve an increase in 



efficiency. To emphasize this periodicity, W will be 



rewritten as Vi in carrying cut the general de velopeme rt . 



nk 

In Eg. £1-9], is a complex exponential, thus 

an examination of ig,£1-9] shews that, in the case when x (n) 
is a complex sequence, a complete direct evaluation cf an 
N-pcint EF‘ I requires (N-1) 2 complex multiplications and 
N (N-1) complex additions. 1 Fig. [9] shews the comparison of 
multiplications reguired by the direct calculation cf the 
EF1 and for the base 2 FFT algorithm. 



4 In most literature, N is assumed to be large, 
sc that (N-1) can he approximated by N . 
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NUMBEH OF MULTIPLICATIONS (xIOUO) 




N (number of sample points) 



Figure 9 



Multiplications Required to 



Compute 



the CFT. 
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t. tecimation-In-Time (DIT) Algorithm 



Ibe general principle behind the FFT is to break 



the original N-pcint sequence into two shorter sequences, 
the DFT's c£ which can be combined to give the DFT of the 
original N-pcint seguence. 

Assuming N is a power of 2, define two 

(N/2) -point sequences x (n) and x (n) as the even and odd 

1 2 

members cf the sequence x(n) respectively, i.e.. 



The N-point DFT of the sequence x(n) can be written as 



x^(n) = x(2n) 
x 2 (n) = x(2n+l) 



N 

n = 0,1, - 1 

M 

n = 0,1, ...,§ - 1 



( 1 - 12 ) 



N-l 



N-l 




n=0 
n odd 



(1-13) 



n=0 



n even 




N/2-1 



N/2-1 



Then using the following identity 



w 2 = (" e j (2-n/N)l 2 = a j[2fl/(N/2)] 




(1-15) 




can written in the form 




N/2-1 



N/2-1 



(1-16) 



uo 



( 1 - 1 ?) 



X(k) = x 1 (k) + w^x 2 (k) 

Each of the sums in Eg. [1-16] is an N/2-pcint 
EFT, the first sum being the N/2 point DFT of the 

even-numbered points of x (n) and the second being the N/2 
point Efl of the cdd-numbered points of the original 

seguence. although the index k ranges over N values, 
k=C , 1 , . . . , N- 1 , each of the sums need only be computed for k 

between C arc N/2-1, since X (k) and X (k) are each periodic 

1 2 

in k with period N/2. After the two DFT • s corresponding to 
the two suns in Eg. [1-16] are computed, they are then 

combined to jield the N-pcint DFT, X (k) as given by 

Eg. [1-17]. 



4 1 



2 . T Jb e Ciitj^Z^Transform 1CZ%1 

The chir p-z- transform (CZT) is an algorithm for 
evaluating the z—transform of a finite duration sequence 
alcng certain general contours in the z-plane. The 
following derivation is based upon a repcrt by Ratiner, 
Schafer, and Eader [18]. 



a. The Derivation of the CZT Algorithm 



In general, the z-transfcrm of a sequence of 
samples jx(n)} is defined as 



00 




n=- co 



Assumminc that the right side of Eg. [1-18] converges for 
seme values cf z, and restricting the evaluation of the 
z-transfcrm to a finite duration sequence of 
sanples, I g .[ 1- 18 ] . 



N-1 




n=0 



The special case cf the z-transform which has 
received ccrsideratle attention is the set cf points ecually 
spaced around the unit circle of the z-plane r i.e.. 
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2 = exp ( j2TTk/N) 

k 



< 1 - 20 ) 



, k=0 



which gives 



X < 2 ) 



N-1 
x (n) e 
n = 0 



- j (2TT/N) nk 



, k=0 , 1 r . . . , N- 1 (1-21) 



By comparing Fg.[ 1-2.1] and Eq. [1-3], one can see that it is 
the discrete Fourier transform (DFT) . The D FT can be 
evaluated acre efficiently by using the FFT rather that the 
CZ1, but the latter allows cne to evaluate Eg. l -19] along 
the mere general ccntcur 



Z 

k 



-k 

= AW 



, k=Q , 1 , ,a-1 



0 - 22 ) 



where a is at arbitrary integer (not necessarily equal tc N) 
and A and W are arbitrary complex numbers defined by 

A = A exp ( j2 tt© ) (1-23) 

c o 



ii = fi exp (.j2tt0 ) 
o o 



(1-24) 



Fig. [10] shews the significance of A^, 9^, 

the z-plane. 



and 0 in 
o 
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Z -PLANE 




figure 10 - Contour of the CZT in the Z-Plane. 
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Icr the special case of the DM, A=1 , M=S, and 
S=exp (— j 2rr/ N ) ; thus the DPT coefficients of a finite 
duration sequence are the values of the z-transform of that 
sane seguence at N evenly spaced points around the unit 
circle . 



The parameter W determines the rate at which 
o 

the contour spirals: if S is greater than unity, the 



contour spirals toward the origin as K increases, and if a 

o 

is less that unity, the contour spirals outward as K 



increases. The parameters A and 2tt9 are the location in 

o o 

radius and angle, respectively, cf the first sample, i.e.. 



fcr k =0. The remaining samples are located alcng the spiral 

contour with an angular spacing of 2^0 . The egui^alent 

c 

s-plane contour of Fig. [10] is shown in Fig. [11]. 
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lfc€ equivalent s-plane contour teginning point 
may be fcund by ttaking the substitution 

sT 

Z = € 



then 



z = A = exp (s T) 
c o 



ci 



s = (1/T) InA 

c 



= a + jw = (1/T) (InA + j2u9) (1-25) 

CO 0 0 



For a genera] point cn the s-plane contour 



-k 

z = exp (s 1) - AW 
k k 



then 



-k 

s ^ (1/T) In (AW ) 

k 



anc finally 

s = (1/T) (lnA-klnW) = s -(k/T)lnW 

k c 

k = 0 , 1 ,M-1 (1-26) 

Thus the spiraling contours in the z-plane correspond to 
straight lines in the s-plane. The rate cf spiraling in the 



46 



z-flane cetemin es the slope cf the lire in the s-plane. 

Ihe task now at hand is to compute Eg. [1-19] 

along a general contour given by Eg. [1-22]. The evaluation 

of X(z) at z=z , will be written using the notation X , 
* k 



N-1 



f = Vx ( 



X. = \x(n)z 

k L k 

n=0 



0-27) 



N-1 

^ 1 -n nk 

X^ = \ x (n> A H , k=0, 1 , . . . , K- 1 (1-28) 



n=0 



where N is the length of the sequence x(n). 

E a king the following substitution, which was 
crigninatea fcy Eluestein [5]. 



nk = nz 4 k.2 - (k-n) 2 

2 



(1-29) 



gives 
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X 

k 



N- 1 




^[n 2 +k 2 - (k-n) ]/2 



n=0 



N-1 



-n n 2 /2 k 2 /2 - (k-n) 2 /2 

[a(n)A W ]W W 



n-0 



(1-30) 



Although Eg. [1-30] appears to be mere complicated than 
Eg. [1-28], it allows efficient hardware implementation. 
Eg. [1-30] it a 5 be viewed as basically a three step process: 



1. The first step consists of a complex 
pre-multiplication forming a new seguence 



y(t) = x (n> A V / ' i ,n=0 (1-31) 



2. The seguence y(n) is then ccnvclved with the 
sequence v(n) defined as 



-n 2 /2 

v ( r ) = W (1-32) 



to give the seguence g (k) defined as 



N-1 

g(k) = y(n)v(k-n) , k-0 , 1 , . . . , fl- h 1 (1-33) 



n = 0 



k 2 /2 

3. The seguence g(k) is then post-multiplied by W 
to cive 
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N-1 



X 

k 



k 2 /2 

w 



y(n)v{k-n) ,k=0,1 






n=0 



(1-34) 



A pictorial repr esen-taticn of the processing operations of 
the CZ1 algorithm is shewn in Fig. [12]. 




A'V 1 / 2 



Figure 12 - Processing Operations of the CZT Algorithm. 
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The advantages cf the CZT algorithm ever the standard 
F5T are summarized below: 



1. N, the number of points in the input sequence, 
doesn’t have to be egual to M, the number of 
points at which the z-transform is evaluated. 



2. Neither N nor M has tc be composite - both can 
be prime numbers. 



3. The starting point in the z-plane and the 
angular spacing of tne z *s are arbitrary. Thus 

k 

the freguency resolution and the range of 
frequencies are arbitrary. 



4. The contour does not have to be a circle in 
the z-plane. The CZT algorithm can evaluate the 
z-transferm along spiral contours inside anc 
outside the unit circle of the z-plane. Ey 
evaluating the z-transform along a selected 
contour, the transfer function of a linear system 
can be adjusted to sharpen or flatten the 
freguency response. The choice of contour is 
depended on which of the system’s poles should be 
emphasized. 



5. If £ = 1, = N , W= exp (-i 2tt/N) , then the CZT can 
be used tc evaluate the DFT, even when N is prime. 



6. fiith the CZT algorithm* interpolation between 
time samples cf a bandlimited functicn can be 
obtained with greater ease with the DFT , and 
inter pclaticn at arbitrary sampling intervals can 
be perfermed efficiently. 
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3. It .§ i i iae Transform 



a. The Prime Transform Derivation 



The fcllpwing derivation is based on a paper by 
Khitehouse, beans, and Speiser [26]. A brief review cf the 
number tbeorj properties used in this derivation is gi?en in 
Appendix B. 



Id the specific case when N is an odd prime. 
Eg. [1-3] can be written as 



N-1 




n = 0 



(1-35) 



to give the (D-C) freguency component and 



X OO-X (0) 



N-1 




j2iTnk/N 



n=1 



,k=1 



, N- 1 (1-36) 



to give the remaining components. 

Since only non-zero values of n and k occur in 
the right hand side of Eg. £1-36], it is possible to replace 
the product tk by a primitive root raised to a sum [19]. 

Proceeding , in the summation cf Eg. [1-36], the 
terms are permuted and the order of the equation changed by 
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th€ following transformations: 



(n)p — ((r^P)) = r n modulo N 

, v , , (1-37) 

(k)p — ((r (k ) p )) = r k modulo N 

n,k=l, a . . , N-l 

where "r" is a primitive root and "p'» implies permuted. 
N-1 0 

Noting that r modulo N = r modulo N f Eg. [1-36] is 

written as 



((r 



OOp\ 



- x(o) 



N-l 



’ I x<n) 



n=0 



((r 



(n)pN 



ex P (-j|2r ((n ^(^P); 



(1-38) 



An examination of Eq, 1-38J, shows that the 

le 

sequence J exp - j ( 2 ti/N) r 1 , 



X - x(0)> is the circular correlation 

((r (k) P)) J 



sequence 
of the 



In matrix notation. Eg. [1-36] may be written as 



i i i 

X - x (0) D = F x (n) (1-39) 



where 




(1-40) 
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f t 

aca X and x (n) are column vectors of size (N-1) derived 
from X ( k ) anc x (a) by deleting, respectively, X (0) and x (0) . 

i 

Ihe matrix F can be factored into three 

matrices 



1 t 

F = P CF (1-41) 

where E is an (N-1) permutaticn matrix, C is an 
[ (N-1) X (N-1) ] circulant matrix, and ? t is the transpose of 
P. 



Ihe elements of the C matrix are 



i 

C = F , n ,k= 1 

r,k (n)p,(k)p 



N-1 



(1-42) 



where "r" is a primitive root of N and '• p *• implies permuted. 
The elements cf the permutation matrix, P, are 



E 

c,k 



s 



(n) p,k 



n ,k= 1 , . . . , N-1 



(1-43) 



where ^ 



J ( d ) E » k 
circulant matrix, C 



trarsversal Hirer 



is the Kronecher delta function. The 
, can be implemented with a recirculating 
whose tap weights are determined by 



h (r) 



w 



(n) p 



, n = 1 , . . . ,N- 1 



(1-44) 
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Ifcus the prime transform algorithm offers more 
sin p lif icaticr than the CZT algorithm because in hardware 
implementation it is possible to eliminate the pre- and 
pest-multipliers . The processing operations maj be 
summarized as 

1. permutation of the input data 

2. circular correlation 

3. permutation of the output Fourier coefficients 
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F. 



hAfiES^BE A BE SCFTHAB2 IMPLEMENTATION 



in tcth hardware and software implementation of theses 
algorithms, the goal is to have the ability to do real-time 
signal processing for the frequency range of interest. In 
spectral analysis applications, signal prosessing takes 
place in real time when the spectrum level cf each freguency 
component is updated at a rate equal tc at least the 
taDdwidth 0 ( B z ) cf the system. 

Depending upon the system, the signal for processing may 
he in any cf three forms, i.e., 

1. Analog-* both time and signal amplitude are 
ccrtinuous 

2. Sanpled analog-* time is discrete and signal 
amplitude is continuous. 

3. Eicital-* time and signal amplitude are tcth 
discrete 



1 • jn alee Spectral Analysis 



a. Scanning Analyzers 



The conventional scanning-f ilter analyzer 
ccnsists cf a single bandpass filter, whose center freguency 
is "shifted" by using frequency- translation techniques. 
Typically a single bandpass filter is employed and the 
filter analyzes cne frequency at a time. 
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Since a single filter is employed and is applied 
successively tc each frequency, it must dwell for a period 
of 1/p at each frequency. Thus (1/p ) X No. of steps equals 
the total analyzing time n T w • The narrower the filter the 
longer is the tctal sweep time. Since each freguency 
component is updated only once per time T, real time signal 
prccessirg is ret possible. 



t. tultifilter Spectrum Analyzers 



Constant bandwidth multifilter analyzers were 
developed for real time analysis. In this type of analyzer, 
the signal passes through many parallel filters 
simultan ecusly , typically 500 contiguous magneto strictive 
filters per analysis range. The filters and the associated 
circuits must possess uniform and constant gain so that a 
fixed aaplitude sine wave of any frequency in the analysis 
range results in a spectrum output which is constant, 
hewever, frequency flatness of better than ±3dB is difficult 
tc maintain, thus the ma gnetcstricti ve multifilter analyzers 
are net used extensively. 



2. Ii£ i t al Spectral Analysis 



The all digital spectrum analyzers employ a computer 
or ccupu ter-like circuitry which is programmed to implement 
the TFT algorithm. Some of the structural factors for the 
design cf such a system are discussed below. 
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TABLE-III 

Comparison of Major Logic Families 
RTL DTL TTL ECL P-MOS CMOS 
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Versatility good fair very good good 



a. Technology 



Easic speed and functional blocks available are 
dictated by technology. Table-Ill shows a comparison of the 
major logic families. Because of speed-complexity 
tradeoffs, tie choice of logic family is usually cot a 
clear-cut decesion. This complexity encompasses many areas, 
but important are available MSI and LSI functional modules, 
system interconnections, multiple function modules, and type 
of multiplier. To gain speed in the FFT algorithm a fast 
multiplier is required in the butterflies, thus the type of 
multiplier is an important consideration. Table-IV is a 
comparison of multipliers. 



TABLE-IV 

MULTIPLIER COMPARISON TABLE 

Computation Time I 





Type of 


1 Component 




(nsec) 




Package Cow 


it 

9x9 


Configuration 


Basic Package 


1 Technology 


16 x 12 


5\ 

OO 


9x9 


16 x 12 


116x8 


Fast trapezoid 


Custom two- 


ECL 


32 


28 


22 


112 


82 


45 


(Fig. 8.35) 


bit adder 
package 










160 


120 


85 


Tree 


Commercial 


ECL 


62 


48 


40 




four-bit 

adder 

package 








85 


160 




85 


Tree 


Commercial 


TTL 


155 


100 


120 




four-bit 
adder 
package 
(2 x 4)-bit 
(commercial) 
array pack- 


TTL 


261 


195 


135 


80 


58 


30 








Two bits at a 


age 

Custom two- 


ECL 


130 


120 


100 


40 


36 


25 


time add- 


bit adder 
















shift 


package 








116 


40 


36 


25 


Two bits at a 


Commercial 


ECL 


210 


140 


time add- 


two-bit 
















shift 


adder 
















Two bits at a 


package 

Commercial 


TTL 


450 


300 


250 


40 


36 


25 


time add- 


two-bit 
















shift 


adder 

package 

















Taken from Ref. [ 17] . 
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t. algorithm 



lh€ structure of the FFT algorithm is a major 
factor in determining whether or not special-purpose 
hardware or a general-purpose computer is used in a 
particular application. Table-V lists the arithmetic 
operations required for different radix algorithms. When 
inplemented with software, the radix-8 algorithm is near 
optimum. However, it lacks the flexability required of a 
general- gui pcse ccnputer and thus is only used in special 
purpose applications. Although the radix-2 requires 
additional computations, its simplicity has made it poplular 
for use in general purpose computers were real time analysis 
is ret a requirement. 



c. hardware Architecture 



When implemented with hardware, the radix-4 
algorithm has received the most attention. Hardware 
in p le men ta tier cf the FFT algorithm has gained speed by 
trading eff flexability. The speed is gained by a complex 
nix cf several types of parallism. One such type is 
pipelining, i.e., the process is broken up' into several 
sequential tasks and execution proceeds assembly— line style. 
Ey pipelining many processing steps can be progress at any 
given time and thus performing real time analysis. Given 
that an N-pcint FFT requires (N/2) log^N butterflies, there 

are four basic ways to organize the processor. 
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* 



Required Arithmetic Operations 
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1 . 



Sequential-* all (N/2) log N butterflies are 

2 

peifcrmed sequentially. 



2. Cascade-* if log^N arithmetic units (AO) are 



available, the algotihm can be pipelined sc 

butterflies 
tie pipe. 



that each AO computes (N/2) 
before advancing data through the 



3. Parallel-Iterative-* if (N/2) AU»s are 
available, the algorithm can be process in a 
leg _ N stage pipeline. This requires leg N 
^ 2 
butterfly times. 



4. Array-* if (N/2)log^N AU»s are available, the 

entire array can oipelined into one execution 
tine . 



Tatle-VI ccnpares these different organizations fer an 
N=1024 F FT ate an execution time of 1 jj s The total execution 
time does net take into account the I/O overhead time. 



TABLE -VI 

Comparison of Processor Organizations 



Machine 

Organization 


Arith. Units 
(Butterflies) 


Total 
Execution 
Time (ps) 


Sequential 


1 . 


5120 . 


Cascade 


10. 


512. 


Parallel 

Iterative 


512 . 


10. 


Array 


5120. 


1 . 



In add ticn tc pipelining, among the types cf parallism are 



1. line overlap cf control and memory functions. 

2. The addition of a small amount of high speed 
memory. 

3. Eicher radix algorithms. 

Tatle-VIl compares some important parameters of a few 
modern signal processing systems. 



3. Sampled Analog Spectral Analysis 

Like the F FT algorithm the CZT algorithm can compute 
the frequency components of a signal mete efficiently than 
hy direct evaluation of the DFT . However, for relatively 
large values of N , the CZT algorithn is less efficient when 
compared tc the FFT algorithm. For this reason interest tc 
digitally implement the CZT algorithm waned while great 
strides were being taken to develcpe FFT software and 
hardware. 



The CZT does offer many signal processing 
advantages, let only if it can be implemented using sampled 
analog techniques. With the developement of Surface 
Acoustic wave (SAW) and Charge Transfer Devices (CTD) , 
interest in the CZT has been re-kindled. The prime 
ad vadvantace cf these devices ever digtal equipment is that 
when the continuous analog signal is sampled the analog 
value is retained, while in a digital system the analog 
value must be quantized which is a source cf error. 

CIE's, which include the Charge Coupled Device (CCD) 
aca Eucket-Erigade Device (2BD) , sample the optical or 
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electrical input and convert it to a packet cf charge whose 
magnitude is prcportinal tc the input. The charge packet 
moves down the CTD stages, shift register fashion, under the 
control cf clock voltages. The charge packets are 
transferred frcn one stage to the next until they readh an 
output stace where the magnitude of the charge packet is 
sensed anc converted to an electrical anlacg signal. The 
deck fregcency and the number of stages transferred 
determine the time delay introduced. 

In c SAW device a mechanical ripple, cr wave, 
travels alone a solid polished surface in much the same 
manner as an ocean wave travels through water. Aluminum 
electrodes are usually deposited on a piezoelectric material 
such as SI quartz, and they act to launch and receive the 
accustic waves. These transducers convert energy between 
electrical ard mechanical domains. 

In a CCD the delay line is tapped and weighted by 
using a split-gate weighting technique which is incorporated 
intc the mask at the time of fabrication or by weighting the 
taps with external resistors. The tap weight in a SAW 
device is phot c-lithographically determined and is 
prcpcrticnal to the length of the interdigitaticn which 
accepts the accustic wave. The pre-and post-multi plies are 
accomplished using wideband four analog multipliers with 
weights generated hy programmable read only memories (EECM) 
and coupled through D/A's or by multiplying D/A's with FROM 
weight stcrace. 

Ihe transversal filter implementation cf the CZT 
algorithm achieves its high speed not by reducing the number 
cf multiplies like the FFT , but by putting the equivalent of 
N parallel multipliers all on the same chip. This reduces 
the number cf computational steps from Nlog^N for the fFT tc 
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U fcr the CZT. Also the N additons per component are 
sin ultanecus . Currently, the CCD CZT can operate in real 
tine at speeds up to 5MHz and the EBD CZT is limited tc a 
sample rate cf a few hundred KHz* For higher sample rates, 
up tc several MHz, SAW devices are attractive. In addition, 
the sampled analog devices reguire lew power and are light 
weight. The advantages of the digital FFT are its 
flexability and accuracy. Programability is being developed 
for the CCD and SAW filters with up tc 128 digitally 
switched taps have been built. 

As rcted earlier, the prime transform offers 
additional simplicity by eliirininating the multipliers 
required in the CZT implementation. The multipliers are 
replaced with analog permuter memories. The commercially 
available serial access memory stores analcg samples as 
charges in an array of MOS capacitors under the control of 
read in and read out shift registers. Currently the maximum 
sample rate fcr this type of device is 5 MHz. 
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Taken from Ref. 1 



II. PAST FOURIER TBANSFORM ALGCBITHMS 



A. KATHEK A1ICAL ILLUSTRATION OF THE DEC IM A1ICN-IN-TIM E 
(CIT) ALGCF HEK 



To illustrate the FFT algorithm, the number cf sample 

y 

points of x (n) is assumed to be a power of 2, i.e., N = 2 . 

In this exaiifle N= 4=2 2 andY=2. 

For S= 4, Eg. [1-3] can be written as 



X(0) = x ( 0 ) W° + x ( 1 ) W° + x ( 2 ) W° + x(3)W° 
X(l) = x(0 )W° + xdlW 1 + x(2)W 2 + x(3)W 3 
X( 2 ) = x(0)W° + x(l)W 2 + x(2)W 4 + x(3)W 6 
X(3) = x ( 0 ) W° + x(l)W 3 + x(2)W 6 + x(3)W 9 

or witten in matrix form 



X(o) 




1 

O 

o 

o 




x(o) 


X(1) 




w° w 1 w 2 w- 3 




x(l) 


X(2) 




w° w 2 w 6 




x(2) 


X(3) 




W° W" 5 w 6 w 9 




x(3) 



(2-1) 



(2-2) 
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Ibe first step in developing the FFT algorithm for this 
example is to rewrite Eg. [2-2] using the following relation 
flea Dumber th€cry[23]. 



lY 

W 



rk 

W 



mcd(N) 



(2-3) 



where [nk icc (N) ] is the remainder upon division of nk by N. 
Hence if = 4 , r = 3 , k=3 



S 1 
fc = H 



(2-4) 



sirce 



nk S 

H = W = exp[ ( — j 2 it/ 4 ) (9) ] = exp (- j 9 tt / 2 > 
= €xp(-jTt/2) = exp[ (- j2 t/4) ( 1) ] 



1 rk mod (N) 

= w = w 

Thus the new natrix is written as 



X(0) 




1111 




x(o) 


X(l) 




1 w 1 w 2 w 3 




x(l) 


X(2) 




2 0 2 
1 W U 




x(2) 


X(3) 




1 w 3 w 2 w 1 




x(3) 



(2-5) 



(2-6) 



The seccrd step is the factorization cf the square 
matrix in fg.[2-€],. As an aid in showing hew Eg. [2-6] is 
factored, a rew notation is introduced at this time, where 
the integers d and k are now represented as binary numbers. 
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n = C , 1 ,2,3 n(n , n ) 
1 0 



CO, 01 



10 



11 



k = C , 1 , 2 ,3 k(k ,k ) = CO, 01, 10, 11 

1 0 

Then writing the tase 10 equivalents 

n = 2n + n ,k=2k + k (2-7) 

10 10 



£g.[1-9] can new t€ written as 
1 1 

X(kl,k 0 ) = N 7x(n 1> n Q )W 



Un^r^) (2^+^) 



(2-8) 



n 0 =0 n ^=0 



Tte important pemt nere, is that the single summation of 
Eg. [1-9] rncst he replaced by " " summations in order to 
enumerate all the hits of the binary representation of n. 



Cc ntinui rg 



(2^+^) (2k 1 +k Q ; 



4n k ' 

IJ - 1 ■ L 



= w 



2 Vl 



W 



w 



( 2k, +k ) 2n n ( 2k, +k ) n 
W 1 u 1 W 1 u 

2 Vl w (2 W n 0 

(2k 1+ k 0 )n 0 



0 



(2-9) 



The term in brackets is equal to unity since 



4n, k. 



1^1 



n i k i 



- j2 4/V 



n l k i 



n i k i 



= l 



( 2 - 10 ) 
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Eg. [2-8] car now te written as 



x < k r k o> = ^ 



V° 



5 



2k n 

x(n lt n 0 )W 0 1 



(2k 1+ k 0 )n 0 



( 2 - 11 ) 



Each of the summations are treated individually and lead 
directly tc the correct factorization cf Eq.[2-6] 
Rewriting the tracketed summation 



'T- 2k o n i 

x 1 (k 0 ,n Q ) = xCr^.riQ) W 

n^O 

Directly evaluating this summation gives 

x 1 (0,0) = x(0,0) + x ( 1 , 0 ) W° 
x 1 (0,l) = x(0,l) + x ( 1 , 1 ) W° 

x 1 (l,0) = x(0 ,0) + x(l,0)W 2 
x 1 (l,l) = x(0,l) + x ( 1 , 1 ) W 2 
re-writing in matrix notation 



x^O ,0) 
x 1 (0 ,1) 
x 1 (l,0) 

_ x i ( 1 , 1 )_ 



10W° 
0 10 
low 2 
0 10 





x(0 ,0) 


) 


x(0 ,1) 




x(l,0) 


l 


x(l,l)_ 



(2-12) 



(2-13) 



(2-14) 
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Now writing the outer summation as 



X 2 ( W = 



( 2k 1 +k 0 ) n 



n 0 =° 



Directly evaluating the summation gives 

x 2 (0,0) = x 1 (0,0) + x 1 (0,l)W° 

x 2 (0,l) = x 1 ( 0 , 0 ) + x 1 ( 0 ,l)W 2 

x 2 (l,0) = x 1 (l,0) + x 1 (l,l)W 1 

x 2 (l,l) = x 1 (l,0) + x 1 (l,l)W 3 

Then in matrix notation 



X 1 

IV) 

o 

O 

— 1 


I 


x 2 (0,l) 




x 2 (l,0) 


( 


x 2 (l,l) 


( 



W° 0 
W 2 0 



1 

1 





x 1 (0,0) 




x 1 (0,l) 




x 1 (l,0) 




x 1 (l,l) 



Ccitininc Egs.[2-12] and [2-15] results in 



X ‘VV = VV V 



or fcritirg it latrix form 



X(0) 




X(2) 




X(l) 




X(3) 





W° 0 
W 2 0 



w 3 



0 

1 

0 

1 



W° 0 



W 2 0 





x(oT 


) 


x(l) 




x(2) 


? 


_x(3) 



(2-15) 



( 2 - 16 ) 



( 2 - 1 ?) 



(2-18) 



(2-19) 
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A signal flew craph of £q.[2-19] is shewn in Fig. [13]. 



X(0) 

X(2) 

X(l) 

X(3) 

Figure 13 - FFT Signal Flcwgraph, N=4 . 

Note that the final results are in bit reversed order with 

respect tc the disered values X(k , k ). This is simply the 

1 0 

scrambli re which results from the FFT algorithm. The 
results ate easily unscrambled by bit reversing, i.e.. 



COMPUTATION ARRAYS 



Data Array Array 1 Array 2 

* 0 (k) Xl (k) x-<k) 




X (CC) 

X ( 1 C ) flips or reverses 
X (Cl) 

x ( 1 1 ) 



X (00) 

to X(01) (2-20) 

X (10) 

X (11) 



The algorithm just described is known as the 
Cccley-Tukey Agcrithm or decimation-in-time (DIT) algorithm, 
since at each stage of the process the input sequence (i.e., 
time sequence) is divided into smaller sequences for 
processing . 

The basic building block of the algorithm is the 
"butterfly" as shewn in Fig. [14]. 
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Figure 14 - Decimatio c-In-Time Butterfly. 

Each butterfly requires one complex multiply (four real 
multipliers) and two complex adds. The importance cf the 
butterfly structure to the FFT flew graph is that only one 
additional nemory location is required to transform an 
N-pcint seguence stored in memory. Thus the intermediate 
results cf the FFT are stored in the same locations in which 
the original cata teing transformed was stored. Since this 
algorithm stores both the input and output sequences in the 
same location, it is called an in-place algorithm. 

As feinted out earlier the results of the DIT algorithm 
are tit reversed. This means that to get a normal ordered 
output, the input data must first be shuffled. Figs. [15a] 
arc [15b] present the canonic signal flew graphs for the DIT 
algorithm for the case N=8. Note that the data is naturally 
ordered ir Fig. [15a] and is in inverse order in Fig. [15b]. 
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figure 15 - FFT Canonic Flew Graphs. 
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B. EECIfilTICS-IN-IBEQUENCY (DIP) ALGOEITHM 



Another popular FPT algcrithm is the Sande-Tukey 
Algorithm or decimaticn-in-f r eguency (DIP) algorithm. It is 
very similar to the DIT algorithm but differs in the 
following ways. 



1. Ihe input sequence fx (n) } is partitioned into 
two sequences each of length (N/2) samples in 
the following manner. The first seguenc 

(x^(n)} consists of the first (N/2) points of 

{ x ( r ) } t where as the second sequence [x (n) } 

consists of the last (N/2) points of (x(n)}. 



2. For a given case in the DIT algcrithm the 
incut is hit-reversed while the output is ir 
natural order, whereas the reverse is true for 
the E IP algorithm. In general, however. both 
the DIT ana DIF can go from normal to snufflec 
data and ?ice versa. 



3. Ihe DIF butterfly is slightly different as 
shewn in Fig.[l6]. The difference being that 
the complex multiplication takes place after 
the add-suttract operation in the DIE 
algcrithm . 




Figure 16 - Decima ticn-In-Fr eguency Butterfly. 
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Taking a lock at the similarities, both algorithms 
reguire cr the order of (N log^N) operations to compute the 

DF 1 . Ecth algorithms can be done in-place and both need to 

perform tit reversal at some place during the computation. 

Figs. [15c] and [15d] present the canonic signal flow 
graphs for the DIF algorithm for the case N = 8. The data is 
naturally ordered in Fig. [15c] and is in inverse order in 
?ig.«15d]. Ihe two most effective methods are those 

illustrated in Figs. [15b] and [15c] since they provide the 
powers c f H ir. the correct order needed for computation. 
This eliminates the need for storage tables. 

The IFT algorithms examined above by no means eshaust 
all the r cssibilit ies. Numerous algorithms have teen 

developed, all with one basic purpose in mind, and that is 
to reduce the computation time. 
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c. 



EXAMPLES C SING THE FET ALGORITHM 



As a test cf the FFT algorithm to compute the DFT, CASES 
1,2, and 3, as described in Appendix A, were used as inputs 
for N = 32. For the LIT algorithm used, N is reguired to be 
highly composite and a power cf 2 • The value of 32 was used 
so that the results of this algorithm could be compared to 
the results cf the CZT and prime transform algorithms. In 
the latter two algorithms N=31. All figures are drawn with 
normalized frequency, but for the sake of discussion, assume 
that the input sinusoids have freguency in Hertz anc the 
rectangular data window "T" is egual to 1 second. 

The fundamental freguency is 

E =1/1=1 Hz 
c 

which is the minimum recognizable freguency compcnert or 
rescluticc. the analyzing bandwidth is 

E = (N/2 ) (E ) = 16 HZ 
max o 

The output is a one-sided spectra consisting of [*N/2) + 1] 
freguency components. Figs. [17] and [18] are the results 
for CASE 1 and those for CASE 2 are shown in Fig. [IS] and 
[ 2C ] . In the results cf CASE 3, the 22 Hz component aliases 
as a 10 E i component. 
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Figure 17 - Magnitude Plot Via FFT Algorithm, Case 
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Figure 18 - Phase Plot Via FFT Algorithm, Case 
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Figure 19 - Magnitude Plot Via FFT Algorithm, Case 
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Figure 20 - Phase Plot Via FFT Algorithm, Case 
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Figure 21 - Magnitude Plot Via FFT Algorithm, Case 
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Figure 22 - Phase Plot Via FFT Algorithm, Case 



III. CH ISP-Z-TR AN SI CR M COrtPOTER SIMULATION 



A. CZT PBOGEAM DE VEICPEUENT 

A computer program was written to simulate the CZT 
algorithm for computing the DFT of an arbitrary waveform. 
This waveform may he real or complex. In carrying cut the 
processing operations, the real and imaginary components are 
represented ty subscripts "r” and "i" respectively. The 
simulation program is a result of the following derivation. 
Proceeding, let 



x-^ input sequence 




—> complex pre-multiplier 



( 3 - 1 ) 



then 



( 3 - 2 ) 




( 3 - 3 ) 



simplify irg 



y = x f - x f 
r r r i i 



( 3 - 4 ) 



y . 

i 



i (X f + x f ) 
r i l r 



( 3 - 5 ) 
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Now letting 



-n 2 /2 



g r y * v 



where "* M denotes convolution 



g = y *v ♦ y *jv + jy *v + jy *jv 
r r ri ir ii 



g = y *v - y *v 
r r r i i 



g = y *jv + jy *v 
i r i i r 



In the program this convolution is simulated 

transversal filter with 2N-1 complex taps v to 

-(N-1) 

where 



-c 2 /2 

v = W , n = - ( N- 1 ) to (N-1) 

n 



and 



« = €Xp (- j2TT/C) 

The third major step in the process, that of 
pcs t- mul tipi ication , is performed fcy letting 



(3-6) 

(3-7) 

(3-8) 

(3-9) 

(3-10) 

as a 
' (N-1) 

(3-11) 

complex 

(3-12) 
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X = gd 



(3-13) 



X=gd + jg d + jg d + j 2 g d (3-14) 

rr r 1 ir ii 

simplifying 

X = g d - g.d. (3-15) 

r r r i i 

X. = j(gd.+g < d) (3-16) 

i r i i r 

Fcr the special case of the DFT, A =1 and 0 =0 in 2g.[l-23], 

o c 

- _y ” n 

K =1 anc 0 =-(l/K) m Eg. [1-24]. Thus A =1, W = exp (- j2 
c c 

rr/N) and Ec.[2-1] # [3-6], and [3-12] can be thought cf as 

complex eifcnential sequences of linearly increasing 
fxeguency. They are represented as follows: 



f = cos (rrn 2 /N) (3-17) 

r 

f = Tsir (Trn 2 /N) (3-18) 

i 

v = cos (TTn 2 /N) (3-19) 

r 

v = sin (ttd 2 /N) (3-20) 

i 

d = ccs (rrk 2 /M) (3-21) 

r 

d * -sic (uk 2 /M) (3-22) 

i 
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Substituting Eg. £3-17] thru [3-22] into Egs.[3-4], [3-5], 
t^“9], [ - - 1 C _ , [j-15], and [3-16], respectively, describes 
the inp lenie ctaticn of the CZT as shown pictorially in 
fig. [23]. 



Y r = x r ccs ( TTn2 / N ) + x.sin(rrnVN) (3-23) 

Y i = j[ -x^sin (Trn 2 /N) + x . cos (Trn 2 /N) ] (3-24) 

c = [y *ccs (rrn 2 /N) ] - [y *sin (rrn 2 /N) ] (3-25) 

r r i 

9. = [ y . * jsin (rrn 2 /N) ] + [jy *cos (Trn 2 /N) ] (3-26) 

11 i 

X = g ccs.^kVil) + g sin(irk 2 /M) (3-27) 

r r i 

X = j[-g sin ♦ g cos (TTk 2 /M) ] (3-28) 

i r i 



Ihe simulation program is listed in Appendix D anc is a 
working icdel of Fig. [23]. 
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Figure 23 - DFT via CZT Algorithm with Parallel Implementation of Complex Arithmetic. 



E. 



EX AM EXES USING THE CZT AlGCfilTHM 



As a test cf the CZT algorithm to compute the DET , CASES 
1,2, and 3 as described in Appendix A, were used as inputs 
for N= 3 1 . The imaginary component of the input is set equal 
to zero since the input consists of only real data. All 
figures are crawn with normalized freguency, but for the 
sake cf discussion, assume that the input sinusoids have 
freguency in Hertz and the data window "T" is equal tc 1 
second. The data window used for all cases is rectangular. 
Thus the fundamental frequency is 



E = _J_ = = 1 = 1 Hz 

o ~ 

(1) N (T ) (31) (1/31) 

s 

which is in effect the minumum recognizable freguency 
component ci resolution. The analyzing bandwidth or 
freguency ranee is 



E = (N/2) (E ) = 15.5 HZ 

max o 



However, since the resolution is 1 Hz, the output components 
consist cf the D-C component, the fundamental, and the first 
14 tanneries fer a frequency range of 15 Hz. 



As in the EFT algorithm, the output of the CZT algorithm 
is a one-sided specta consisting cf [ (N/2) +1 ] integer 
freguency components. Figs. £24] and [25] are the expected 
results fer CASE 1 and Figs. [26] and [27] are the results 
for CASE 2. In the results of CASE 3, it should be noted 
that the aliasing frequency of 22 Hz fclds back as a 9 



88 



Hz, whereas in the FFT it folded hack as a 10 Hz component. 
This occurs tecause N is odd in the CZT where as in the FFT , 
N is required to he even and highly composite. For the FFT 
alccritha N =22. 
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Figure 24- Magnitude Plot Via CZT Algorithm, Case 
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Figure 2 5 - Phase Plot Via CZT Algorithm, Case 
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Figure 26 - Magnitude Plot Via CZT Algorithm, Case 
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Figure 27 - Phase Plot Via CZT Algorithm, Case 
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Figure 28 - Magnitude Plot Via CZT Algorithm, Case 



180 




o 




.00 



VQ 



+ 



— u 

o 

o\ 



o 



o 



(SOOaSep) 0SBLtc£ 



o 

ON 

I 



o 

00 



95 



Figure 29 - Phase Plot Via CZT Algorithm, Case 



IV. £1111 E TRANSPORT COMPUTER SIMULATION 



A. PRIME TEHNSPORM COMPUTER EBOGRAM DEVELCPEMENT 



So as tc make the mathematics clear, a simple example is 
given tc illustrate the required matrices cf the algorithm. 

Consider the case N = 5 and let W = exp(-j2 tt /N) . 
Evaluating the DPT equation directly results in the 
following N eguaticns. 



X(0) = x(0)W° + x(l)W° + x(2)W° + x ( 3 ) W° + x(4)W° 

X(l) = x(0)w° + xCDW 1 + x(2)W 2 + x(3)W 3 + x(4)W 4 

X(2) = x(0)W° + x(l)W 2 + x(2)W 4 + x(3)W 6 + x(4)W 8 (4-1) 

X(3) = x(0)W° + x(l)W 3 + x ( 2 ) W 6 + x(3)W 9 + x(4)W 12 

X(4) = x ( 0 ) W° + x(l)W 4 + x(2)W 8 + x(3)W 12 + x(4)W 16 
Written in matrix form 



X(k') =[pj L x(n)] 



’x(o) 




1 

O 

o 

o 

o 

0 

1 




x(0) 


x(l) 




w° w 3, w 2 w 3 w 4 




x(l) 


X(2) 


- 


w° w 2 w 4 w 8 w 8 




x(2) 


x(3) 




w° w 3 w 8 w 9 w 12 




x(3) 


X(4)_ 


1 


[W 0 w 4 w 8 w 12 w 18 


I 1 


_x(4)_ 
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-he matrix {F} is written using the relationships 
w nk = w nk mod N 
W° = l 

For example W 9 = W 4 since for n=3,k=3 

w nk = w 9 = exp [(.j^ )(9) ] = exp(-ji|S) 

= exp(-j|^) = exp [(-j4?) (4)] = W 4 

therefore 



X(o) 




11111 




x(0) 


X(l) 




1 w 1 w 2 w 4 




x(l) 


X(2) 


= 


l w 2 w 4 w 1 w 3 




x(2) 


X(3) 




1 w 3 w 1 w 4 w 2 




X(3) 


i_X(4) 


| 


1 w 4 w 2 w 1 


1 


_x(4) 



Then deleting X(0) and x(0) , Eq. [4-4 ] is 



X(l) 




l" 




W 1 


w 2 


w 3 


w 4 




x(l) 


X(2) 




1 




w 2 


w 4 


w 1 


w 3 




x(2) 


X(3) 


- x(0) 


1 


“ 


W 3 


w 1 


w 4 


w 2 




x(3) 


_X(4)_ 








w 4 


w 3 


w 2 


w 1 




x(4) 



X* D F* x(n) 

For N=5 and r=2 , there is a one to one mapping of 
integers n,k to the integers (n)p,(k)p, i.e., 

2 1 = 2 mod 5, 2 2 = 4 mod 5, 2 3 = 3 mod 5 , 2 4 = 1 mod 5 



n,k 


1 


2 


3 


4 


(n)p, (k)p 


2 


1 4 


3 


1 



(4-3) 



(4-4) 



(4-5) 



the 



(4-6) 



97 



Matrix {F } is factored into three matrices. 



F = P CP 



The elements of the [C} matrix are found by using the 
relationship given by Eq.[l-42l. Thus 



C = 



W 3 



W J 

w 1 



vr w 

,2 T . T 4 



w 1 - 



vr vr 
w 4 



(4-7) 



vr vr w - 7 w 

The elements of the permutation matrix are given by 
Eq.[l-43]„ 

r 0 1 0 0 

0 0 0 1 

0 0 10 

1 o 0 0 



(4-8) 



The elements of the transpose of P are 



0 

0 

0 

1 



0 1 
0 0 
1 0 
0 0 



(4-9) 



A functional diagram of the prime transform algorithm 
for computing the Fourier coefficients of real data is as 
shown in Fig.i 30 . The switch is in the up position for the 
first data sample and is down for the remaining (N-l) sam- 
ples . 



S8 



sample and hold 
(hold for N time units) 



x 




Cti 

0 rH 

a i 
z 
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Figure 30 - Functional Diagram of the Prime Transform Algorithm. 



B. EXAME IBS GSING THE PRIME TRANSFORM ALGORITHM 



Tc simulate the operations of the prime transform 
algcritha, a computer program was written for the odd prime, 
N=31. A listing of the program is given in Appendix E. 

Once again Cases 1,2, and 3 were employed in order tc 
evaluate the prime transform algorithm. The fundamental 
freguency and frequency range are the same as that cf the 
CZT algorithm for N=31. The results are shewn in Figs. [31] 
thru [36]. As can he seen, the results are the same as 
these computed by the CZT algorithm. 
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Figure 31 - Magnitude Plot Via Prime Transform Algorithm, Case 
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Figure 32 - Phase Plot Via Prime Transform Algorithm, Case 
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Figure 33 - Magnitude Plot Via Prime Transform Algorithm, Case 
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Figure 3^ - Phase Plot Via Prime Transform Algorithm, Case 
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Figure 35 - Magnitude Plot Via Prime Transform Algorithm, Case 
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Figure 36 - Phase Plot Via Prime Transform Algorithm, Case 



V. CZT SENSITIVITY ANALYSIS 



This cheater addresses the sources cf error ir Loth 
digital and sampled analog systems with special attention 
given tc the sampled analog implementation of the CZT 
algcrith s . 



A. DIGITAL JIT SYSTEM 



In a digital FFT system, the sources cf error nay be 
summarized as fellows; 



1. Cuantization of the data at the input A/£ 
converter . 



2. Inaccurate transformation due tc errors in 

k 

finite representation cf the coefficients W 

K 

in the butterflies. 



3. Ecundoff noise in truncating the result of a 
multiplication and scaling the data tc 
prevent overflow. 



An rms quantization error in the A/E converter can be 
defined as 



-b 1 

Q = 2 / 12 

e 



(5-1) 
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where fcl is the nuiber of bits used to represent the input 
cata. This ncise dees not scale with signal size and 
dciinates cnly at low signal levlels. 

The errer due to the inexact representation cf the 
]£ 

multiplier coefficients W can be expressed as the ratio of 

N 

mean sguare output error to mean sguare output signal to 
give 



E 

e 



[ lcg 2 N 

\ r 



2 ” (2) (t3) 



(5-2) 



where b3 is the number of bits used to represent the FFT 
coefficients. This eguation predicts that the variance of 
the error increases at a rather slow rate with N. 



A measure of the rms noise output 
output due to errors produced by 
trrncaticn to prevent overflows in the 
by 



to the rms signal 
roundoff and signal 
butterflies is given 



r is (error) 



f— -b2 

VN2 (0.3) 



rms (signal) 



rms (input) 



(5-3) 



where b2 is the nuiber of bits used to represent the result 
cf the operation. This upper bound of the error- to- signal 
is seen to increase a s~\fW, and is the mest iiportant source 
of error in a digital FFT. 
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E. 



SAMP LIE ANALOG CZT SYSTEM 



Instead cf errors resulting from registers containing a 
finite word length , as in the digital EFT system, the output 
accuracy cf tie sampled analog CZT system is limited by the 
CCE's and SAfc's sensitivity to the environent (temperature, 
humidity, etc.,) and fabrication imperfections. The sources 
cf error cf a CCC CZT system using pre-* and post-multiplying 
coefficients stored in either a ROM or PROM are as follows 

1. Thermal noise-> this source cf error is 
analogous to guantization in a digital FF1 
because it generates an error which is 
independent of signal level. 

2. The errcrs introduced by inaccuracies in the 
pre-and post-multipliers. 

3. The errors introduced by the inaccuracies in 
the tap weights of the transversal filter. 

4. The error due to mismatch and nonlinearities 
in the differential current amplifiers. 

5. The error due to mismatches of the 

p cst-multiplier outputs. 

The ericrs addressed in this study where those due to 
inaccuracies in the pre- and post-multipliers and the 

inaccuracies in the tap weights cf the transversal filter. 
The errors in the pre- and p cst- multipliers were modeled as 
being a percentage "p M of the theoretical coefficient value 
m (n) , i. e . , 

V = p[ m (n) _] (5~4) 
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and to be uniformly distributed on (-V2, */2) . The 

transversal filter was modeled as a tapped analog delay line 
fcith the tap weights determined by the value of the external 
resistors. The error in the value of the tap weight was 
modeled as being a percentage "p" of the theoretical tap 
eright value k (n) , i.e., 

P = (5-5) 

and tc be uniformly distributed cn (-p/2, p/2). These error 
models were included in the simulation program [Appendix D]. 



C. SENSITIVITY TESTS 

Three different input waveforms were used to evaluate 
the perfcrmace of the CZT algorithm for the following 
percentages 

p = Cl, 4%, 10%, 20%, 30%, and 40% 

In all three tests N=512 and a rectangular data window equal 

tc 4 seconds was used. The value of 4 seconds was used for 

ease of plotting the output spectrum. Thus the frequency 

resolution is f/f = 0.25 and the frequency range is f/f = 
s s 

64. C. 



The first waveform consisted of the impulse function 
shewn in fig. [37]. The time impulse responses of the cosine 
and sine filters of the convolver for the different 
percentages are shown in Eigs.[38] thru [50]. For 
p=C%,Fig .[ 5 1 ] shows the expected results, i.e., equal 
magnitude frequency components covering the entire frequency 
ranee. As tte percentage of error increases, examination of 



the figcres ieveals that the average vale cf the magnitude 
of the frequency components is greater than 1.0. 
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Figure 37 - Impulse Input (normalized A/N) . 
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Figure 38 - Time Impulse Response of the Cosine Filter ;N=512 , p=0$. 
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Figure 39 - Time Impulse Response of the Sine Filter jN= 512 , p=0%. 
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Figure 40 - Time Impulse Response of the Cosine Filter ; N=512 , p=4$ 
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Figure 4l - Time Impulse Response of the Sine Filter ;N=512 , p=4$. 
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Figure 42 _ Time Impulse Response of the Cosine Filter ;N=512 , p=10^. 
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Figure ^3 - Time Impulse Response of the Sine Filter ;N=512 , p=10$. 
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Figure 44 - Time Impulse Response of the Cosine Filter j N=512 , p=20$. 
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Figure 45 - Time Impulse Response of the Sine Filter ;N=512 , p=20^. 
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Figure 46 - Time Impulse Response of the Cosine Filter ;N=512 , p=30$* 
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Figure 47 - Time Impulse Response of the Sine Filter ?N=512 , p=30^. 
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Figure 48 - Time Impulse Response of the Cosine Filter ;N=512 , p=40$. 
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Figure 49 - Time Impulse Response of the Sine Filter ;N=512 , p=40%. 
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Figure 50 - Magnitude Plot of Input Impulse, p=0%. 
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Figure 51 - Magnitude Plot of Input Impulse, p= 4 $. 
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Figure 52 - Magnitude Plot of Input Impulse, p=10$. 
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Figure 53 - Magnitude Plot of Input Impulse, p=20$. 
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Figure 5^ - Magnitude Plot of Input Impulse, p=30^. 
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Figure 55 - Magnitude Plot of Input Impulse, p-40 



Fig. £57] was the second waveform used fcr computinc the 
DF1. The magnitude of the one-sided spectra is shewn in 
Figs. £58] thru £63], Degradation of the output spectrum was 
never sc severe that the (sin x) /x function was not 
reccgnizatle , hut as pointed out above, the average value of 
the magnitude shows a moderate increase as the value cf "p" 
increases . 

The third input waveform, shown in Fig. £64], consisted 
of the fcllcwirg components. 





Amp. 


f/f 

s 


Phase 

Shift 


sine 


1.C 


14.125 


0O 


cosine 


-0. 6 


15.0 


Oo 


sine 


1. C 


16.0 


0O 


cosine 


0 . € 


16.25 


0o 


sine 


-0.4 


17.5 


600 



The normalized one-sided spectra and phase plots are 

shewn ir Fig. £65} thru £76], In Fig. £65], the spreading of 

the 14.125 component is due to leakage, the 16.0 and 16.25 

components appear to be indistinguishable, however, this is 

a miss-leading result caused by the continucus curve plot. 

If a point graph had been used, the two components would be 

readily identified, an expected result since the freguency 

res<luticr is f/f - 0.25. 

s 

Examination of the phase plots reveals that the 
degradation is mere severe than that of the magnitude plots. 
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ev€D for snail values of "p". As in the previous exaiiples, 
the magnitude plots also exhibit a general increase in the 
component magnitude except for the f/f = 16.0 component. 

c 

This component shows a very distinguishing decrease as "c" 

increases. Ej re-examining the magnitude plots of the 
impulse waveform, it is obvious that the f/f = 16.0 

c 

component will sho-w a decrease as "p" increases. These 

plcts also show that a more severe degradation cf the 

one-sided spectra would result for a similar set of 

ccnpcnents near f/f = 48.0. 

s 
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Figure $6 - Rectangular Pulse Input. 
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Figure 57 - Magnitude Plot of Rectangular Pulse Input, p=0^. 
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Figure 58 - Magnitude Plot of Rectangular Pulse Input, p=4$. 
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Figure 59 - Magnitude Plot of Rectangular Pulse Input, p=10$. 
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Figure 60 - Magnitude Plot of Rectangular Pulse Input, p=20fi. 
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Figure 6l - Magnitude Plot of Rectangular Pulse Input, p=30$. 
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Figure 62 - Magnitude Plot of Rectangular Pulse Input, p=40$. 
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Figure 63 - Sinusoidal Input. 
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Figure 64 - Magnitude Plot of the Sinusoidal Input, p=0 %. 
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Figure 65 - Phase Plot of the Sinusoidal Input, p=0$. 
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Figure 66 - Magnitude Plot of the Sinusoidal Input, p=4$. 
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Figure 67 - Phase Plot of the Sinusoidal Input, p= 4 $. 
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Figure 68 - Magnitude Plot of the Sinusoidal Input, p=10^. 
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Figure 69 - Phase Plot of the Sinusoidal Input, p=10^. 
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Figure 70 - Magnitude Plot of the Sinusoidal Input, p=20%. 
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Figure 71 - Phase Plot of the Sinusoidal Input, p=20%. 
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Figure 72 - Magnitude Plot of the Sinusoidal Input, p= 30 %. 
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Figure 73 - Phase Plot of the Sinusoidal Input, p=30$. 
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Figure 74 - Magnitude Plot of the Sinusoidal Input, p= 4 o^. 
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Figure 75 - Phase Plot of the Sinusoidal Input, p=40$. 



VI. CO NC LOS 10 NS AND BECCMMENEAIIONS 



The EFT algorithm is a well established method for 
confuting the DF1. The chirp-z-transf orm algorithm can also 
compute the DFT and since the bulk cf the required 
computations can be performed by a transversal filter, it 
lends itself naturally to implementation with sampled analog 
devices. The prime transform algorithm has demonstrated its 
ability to compute the DFT and like the CZT algorithm, can 
be implemented with sampled analog devices. The fact that 
the multipliers reguired by the CZT algorithm are replaced 
by permuting memories in the prime transform algorithm may 
or may net be an advantage. It will depend upon the speed 
and dynamic iange reguired for a particular application. 

This paper evaluated the errors in the weights cf a 
transversal filter for weighting provided externally tc the 
CCE by resistors. Preliminary results indicate that the tap 
weights may vary as much as 20% before the freguency 
components are greatly distorted. However, additional 
analysis needs to be conducted before conclusive results are 
presented. In addition to resistor weighting, a split-gate 
weight technigue has been developed which weights the taps 
or the transversal filter during fabrication. For this type 
cf weighting, quantization errors are introduced. This 
suggests that a different error model will re reguired to 
ccr cuct a sensitivity analysis. 
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APPENDIX A 



ALG0BI1HM TESI DATA 



Several test waveforms were computer generated to 
evaluate the performance of the various algorithms used to 
compute the EET. 



Case 1- 


> This 


waveform 


consisted cf 


the 


following 


ccmpcneEts : 


Amp . 


f/f 

c 


Phase 












Shift 






E-C 


0.8 


0.0 


0 o 






cosine 


1.0 


2.0 


0 ° 






sice 


2.0 


6.0 


0 ° 






cosine 


1.0 


14.0 


450 






The intent 


cf this 


case was 


to see how well 


the 


algorithms 



identified the magnitude and phase information cf each 
frequency component. 
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Case-> This waveform consisted of the following 
ccipcnents 



Amp . 


f/t 

s 


Phase 

Shift 




t-C G . 6 


0.0 


00 




cosine 1.0 


2.0 


0O 




sine 2 . C 


6.0 


0O 




sice -1.5 


10.5 


00 




ccsire l.C 


14.0 


450 




This case illustrate 


s what happens when 


leakage is present. 


The data kindcw truncates the 


f/f = 10. 

s 


5 component at a 


non-integer multiple 


of the 


wavelength. 


which is the cause 


of leakage. 


Case 3-> This 


waveform 


consisted 


of the following 


cca pcnen t s 


Amp . 


t/t 

s 


Phase 

Shift 




i-i 

i 

n 

o 

CD 


O 

O 


00 




cosine l.C 


2.0 


0o 




sine 2. C 


6.0 


0o 




cosine 1.0 


14.0 


450 




ccsire -2-C 


22.0 


0O 
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thus 



This waveform has a component greater than (f )/2 , 

s 

this case was used to illustrate the aliasing phenomenon and 

t fc € effect K (being odd or even) had on the folding cf this 
component into the analyzing bandwidth. 
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APPENDIX B 



EBIEF REVIEW OF INTEGERS AND INDICES 



The following definition and theorems are provided for a 
tetter understanding of the derivation of the prime 
transfers algorithm 



A . DEFI MUCKS 



Unit-* a unit is an integer that divides every integer. 
Since +1 and -1 divide every integer, they are both units. 

Null element-* it is an integer that divides only 
itself. Kc integer different from zero is a null element. 
Zero is the null element of the rational integers. 

Associates of an integer-* they are the results of 
multiplying the integer by the units. 

Prime-* it is an integer, not a unit, that is divisible 
by cnly its associates and the units. This definition 
implies that the greatest common diviscr of a prime "p" and 
an integer "a" is 1, or the positive associate of ,, p n . 

Composite-* it is an integer that is not the null 
element, a unit, or a prime, e.g. the integer 256 is said to 
be highly composite while the integer 563 is moderately 
composite . 
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Relatively prime-> two or more integers are prise to 
each ether cr relatively prise, if their greatest common 
divisor is +1. The integers 6,-9, 14 are relatively prime. 

Tctiert cf m-> denoted p(m) is the number of positive 
integers less than "a" and relatively prime to "m" • 



E. THEOREMS 

1. Iwc integers "a** and ,, b" are said to be congruent 
modulo and integer "a” if "m" is an integer divisor cf a-b. 
This relaticr between ?, a*' and M b" is most often denoted by 

a s b mod <m) (E-1) 

and is read "a is congruent to b mod (m) 11 Eut for the ease 
cf notation, it is denoted as 

( (a) ) = b mod ( m) (E-2) 

in the derivation cf the prime transform algorithm. 

If "a" is congruent to "b* 1 modulo (id) and 0<b<m, ,, b ,, is 
called the residue of a module (m) . If "m" is net an 

integer divisor of a-b, "a” and "b" are said to be distinct 
mcdulo m. 

2. If "a” and "m" (<0) are relatively prime integers, 

then 



jZ(m) 
i (a 



)) 



1 mod (m) 



<E-3) 
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Corollary of 2 . If “p” is a prime which is not a 
divisor of at integer "a", then 

p-1 

((a )) 55 1 mod (p) (E-4) 

If "a' 1 and "i' 1 (>0) are relatively prime integers, and ,, u M is 
the least positive integer for which 

u 

< (a )) = 1 mod (u) (E-5) 

"u" is called the exponent to which "a” belongs modulo (m) . 
If the exponent to which '’a" belongs module <m) is #{m) , "a" 
is defined as a primitive root modulo (m) (or a primitive 
icct of ii) . 

3. The only integers which possess primitve roots are 
n r 

2,4,p ,2p , where "p" is an odd prime. If "m" does have a 
primitive root, it has 0C0(m)) distinct roots module (m) . 

Let "p" be any prime, and let "r" be any primitive root 
cf "p". lo each integer "a" relatively prime to "p" there 
corresponds a unique integer M i" such that 

((a)) = r 1 mod (p) (0<i<p-1) (E-6) 

The integer "i" is called the index to the base "r" cf a 
modulo (p) . This is written as 



i 



ind a 
r 



(E-7) 



Properties of Indices 



(a) ind 1 = 0 
r 


<E-8) 


(t) ind (-1) = <p-1)/2 
r 


(E-9) 


(c) ind (ab) = ind a + ind b mod (p-1) 
r r r 


(B-10) 


n 

(a) ind a = n ind a mod (p-1) 
r r 


(B-11) 


The power residue of an integer "b 1 * to the base 
simply the integer "s" such that 


"r" is 


b 

((s)) = r med (p) (0<s<p) 


(B-12) 


Since the power residue of ind a is congruent 

r 


tc a 



module (p) , the concept of indicies and pewer residues is 
analogous tc that cf logarithms and antilogrithms. 
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APPENDIX C 



PRGGRAM LISTING OF THE F FT ALGORITHM 



ThlS PROGRAM USES THE CEC IMAT ICN-IN-T IMS ( CCCL E Y-TUK E Y ) 
F FT ALGORITHM TC COMPUTE THE DISCRETE FOURIER TRANSFORM 
( C FT ) COEFFICIENTS OF A COMPLEX INPLT. 



C IMENSICN 
C 1MENSICN 
C1MENSICN 
C 1MENSICN 
CINENSICN 
DIMENSION 
CCMPLEX 2 
COMPLEX *8 
INTEGERS 
P EAL*4 



T V ( 10 2 4 ) 

F V ( 5 14 ) 

V ( 1024 ) 

AMPU024) » P HA SE ( 600 ) 
S(1012) tINV(1012) » M ( 3 ) 
TIMEUOC) ,XS(600) 

4(1024,1, 1) 

I TB < 1 2 ) / 12^0/ 

RTB ( 23 )/ 28*0.0/ 



RE AC IN THE POWER OF TWO "KK" WHICH DETERMINES THE NUM- 
EEP CF CATA PCINTS ANC CCMPLEX FREQUENCY COMPONENTS ANC 
THE LENGTH OF THE CATA WINCCW «T" IN SECONDS « 

GEAC(5,6)KK,T 
t FC P M A T ( 15, E10.7) 

N 1=2**KK 
N 2= 1 
N 3= 1 
N h = M/2 
T S « T/M 

CCMPLTE THE SAMPLED VALUES OF AN ARBITRARY REAL WAVEFORM. 
CC 215 j = 1 , N1 

5s(j)’= 0.8 + C0S(6. 28318 * 2. * K * T S ) + ( 2.0 )*S IN < 6. 
i2i z IS * 6.0 * K * TS) + 

$ C^S ( (6.28318 * 14. * K * TS) + 0.7E5ES) 

S- (2.0)<CCS(6.23318 * 22.0 * K * T S ) 

215 CONTINUE 
N (1 1 =KK 
V (2 )=0 
M (2) = 0 

CALCLLATE THE FUNDAMENTAL FRESLENCY ANC ITS FAFMCNICS. 
CALCLLATE THE TIME AXIS VALUES. 

FLNC = (l./T) 

F V ( 1 ) =C .0 

CC 20 J s 2 »NH 

fv( j)=fv(l-ii+fl;nc 

2C CCNTINLE 
TV ( 1 )-0 .0 
CC 20 X s 2 t N 1 
TV <K)=TV(K-1)+TS 
2C CCNTINLE 
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JET ThS C0FFLEX "A" R ATR IX TC ZERC. 

CC 2 11 = 1, N1 

. CC 2 ll = i InI 

CC 3 1 2 = 1 , N2 
CC 3 13 = 1, N3 

TRANSFER THE SAMPLED INPUT OATA TO THE COMPLEX -A" PA- 

3 CtHifoi' ' 

JNFLWATa! DISCRETE FCUPIER TRANSFORM OF THE SAMPLED 

Cm HAPM(A f M, INVfSf-l , I FER P ) 

FCNENTsT E TH£ ,VAG,NITuce PHASE CF THE FRECUENCY CCM- 

CC 33 I 1 = 1 f NH 
Z * A(I 1 , 1 , 1 ) 

AFEAL = REAL(Z) 

ZINAG = AINAG(Z) 

IF(AB$(ZREAL).LT.Q.G01) zreal = c c 

i^Fti i > 2 i"c§esVzT C * 001 ’ ZIM4G - 

FMSEdlJ = 1C0.0 

TcMS ? A. L *i?-2-5-^^ C - Z I yAG » Ge - 0 -0)PHASE<Il) 

r e 7§c^*rc , §*2*^£*£J MAG * LT * 0 * 0 >Pt- a SE< II ) 



IMZRIAL.GE.O.O.ANoIziMALscroIcIprASFlu = 
lf<ZREAL.Li.O.O.AND.ZIMAG.EC.O.(3)FI-ASE(Il) = 
it !£££§?< I1J*NE.1GG.0) GC TC 33 
, . f rsT f V, 1 C = (ATAN21ZIMAG, ZREAL) ) * 57.2S57 
z : LLNTINLE 

CL7PLT THE COMPUTED CATA. 

WFI7E(6,47> T,TS,M,FLNC 

47 FCFMAT(/,* SAMPLE TIME = SFT.Af* SECONDS*//* 
$' SAMPLE INTERVAL = *,F6.4 1 SEfcNCW/t 
i' NUMBER CF SAMPLES = M4// f 
S' FUNDAMENTAL F R ECUENiC Y = «,F6.1,* HZ') 

IT E (3 ) = 8 
176(A) = 5 
176 < 12 ) = 1 
N F = N 1 
PTE (1) = C.l 
P 7 £ ( 2 ) = 2.0 

FLC7 THE SAMPLED INPUT WAVEFORM. 

CALL CRAWP <NH,TV,XS,I7B,RTE) 

176 (2) = 2 
N F = Nl/2 
P7E(1) = 2.0 
F 7 6 (2) = 0.2 

FLC7 THE MAGNITUDE OF THE CNE-SIDEC SPECTRA. 

CALL CFAWP(NH, FV, AMP ,ITB,RT5) 

PTE ( 2 ) = 90.0 

PICT THE PHASE GF THE CNE-SIDEC SPECTRA. 

CALL CRAWP(NH,FV f PHASE, ITe,RTB) 

S TC P 
ENC 



90.0 
-90.0 

c.o 

160.0 
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APPENCIX C 



PROGRAM LISTING CF THE CZT A LGQRI TPP 



*£S^?!=!IS S Cl SC P ET E FPECU5NCY COMPONENTS 

cc^oJ I I?^?X-S 92 PLex INFUT mvefcrp ey tfe chir°-z- 

IrANiFCRM ALGORITHM. 

CCJJCN X(1024) ,Y< 1024), XW( 512), YW< 512), XS< 512) 

CCNNCN >?( 5U ,XP( 512) tYP( 512) ,FREC( 512 ) ,VR <5121 
CCM NCN Wl{512),GR(512),GI(512),AMp(5l2),TIME<512) 
rrttrl^ £ f E < ?1 \ * A 1 1 { 5 1 2 i » F y T , N , P N , S P L , P R P W , T R W » P C P V 
CCMMCN FCTlyPCT2,FCT3,PCT4, FCT5,FCT6 
INTEGERS ITB(12)/12*0/ 

R E A l* 4 RTe(28 )/28*C.O/ 

PE*C THE PARAMETER "II", WHICH CETEPPIN6S IF CATA IS TO 
EE F E £ C IN CP CALCULATED WITHIN THE PROGRAM. IF "II” IS 
ECLAL TC "C", THE PRCGRAM WILL GENERATE A CCMPLEX SINL- 
SCICAL INPUT FOR A SINGLE FREQUENCY • IF ”I1» I< NCT 
EQUAL TD ZERO, AN ARBITRARY REAL INFUT MAY EE F £ AC IN ANC 
THE IMAGINARY PART (YS) IS SET TG ZEPG. 

IFLN = 1 
FEAC(5, 11C) II 
11C FCFMT(I5) 

"N ” IS THE NUMBER OF TIME ANC FREQUENCY SAMPLE PCINTS. 

"F" IS TP? FREQUENCY TC BE USEC IN SUBROUTINE WAVE ANC 

"T" IS THE LENGTH OF THE INPUT CATA WINCCW IN SECCNCS. 

" S F L " CETEFMINES THE PATE CF SPIRALING CF THE CONTOUR 
IN THE Z-PIANE. IF «$PL» IS EQUAL TC ZERO, THE CCNTCLR 
IS A UNIT CIRCLE IN THE Z-PLANE WHICH CCPRESFCNCS TC 
CALCULATING THE DISCRETE FCURI5R TRANSFORM. 
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R E A C ( 5 , 105) F,T,SPL,N 
FCF NAT ( 2610.7 ,1 5 ) 



MAC IN THE PARAMETER "FRRW" WHICH CETERMI NFS WHETHER CR 
NCT FANDCN VALUED FR E-MU LT IPLY INC COEFFICIENTS WILL EE 
USEC IN THE COMPUTATION. PARAMETER "TRW" DETERMINES 
WHETHER CR NCT RANDOM VALUEC TAPPING WEIGHTS IN THE CCN- 
V C LL T I CN FILTERS WILL EE USED IN THE COMPUTATION. 
PARAMETER "PORW" DETERMINES WHETHER CR NCT RANCCN VALUEC 
FCST-MITI FLYING WEIGHTS WILL EE USEC IN THE CONFUTATION. 

F $ A C ( 5 , 101) PRRW , TRW, PCPWI 
1C1 FCFNAT(2E10.7) 

"FCT1" ANC "FCT2" DETERMINE THE RANGE CF VALLES THE 
FPE-M.ITIFIYING COEFFICIENTS CAN CEVIATE FROM A NOMINAL 
VALLE. "PCT3 " AND " FCT4" DETERMINE THE RANGE CF VALLES 
Tpc FILTER TAPPING WEIGHTS CAN CEVIATE FROM A NOMINAL 
VALLF. "P CT 5 " AND " FCT6" DETERMINE THE RANGE CF VALUES 
FCT- MULTIPLYING COEFFICIENTS CAN DEVIATE FR CM A NOMINAL 
VALUE. 
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ICS FCT2 ’ PCT3 ’ PCT <’ PCT5 ’ PCT4 

FFEC(l) = 0.0 
F f 1 = N 

ec IF(U.EC.O) GO TC 5 

.! E£u “ ,c IER0 
CC 40 M = 1,N 
.. = M + 1 
YS(M> = C.C 
F F EC ( J I = M * ( 1./7 ) 

4C CONTINUE 
GC TC 9 

E/NFIEC *N » T T(5ES S ° F ARelTpilRY COMPLEX V» d V E F C R M 
5 CUL WAVE 

9 cm ffm.lt 

cm xslm 
C / L L T 4 F V» 

cm ccnv 

C/Ll FCMJlT 
C/U AMFFA 

CLTFLT TFE COMPUTED CATA. 

WFITE(6»320) IRUN » FCT2 
220 FCFMAT< / ,2X, 'THIS PUN IS NUMBER ',13, 

J* PNC THE PERCENT DEVIATION IS <,F2.2) 

HE (2 I = 6 
IT E (A ) = 5 

it e < 12 ) = i 

FTE (i) “ 0.1 
FIE ( 2 ) * 2.0 
MLMFTS = N 

FLCT TFE SAMPLED INPUT WAVEFORM. 

C/LL C P Atop ( NUMPTS .TIME » XS » I TE »RT E ) 

FIE (2 ) = 0.4 

FLCT TFE IMPULSE RESPONSE CF TFE COSINE FILTER. 

C/LL CRAWFtNLMPTS,TI''E,X,ITR,RTE) 

FLCT TFE IMPULSE RESPONSE CF TFE SINE FILTER. 

C/IL CF AFF (NUMPTS.TIME, V, ITE.RTE ) 

NLMFTS = N/2 
FTE (1) = 2.0 
FTE (2) = 0.2 

FLCT TFE MAGNITUDE CF THE FREOUENCY COMPONENTS. 

C/LL CF/NP(NUMPTS , FPEC.AMP, ITS, PTE ) 

PTE <2 I = 90. C 

FLCT TFE PFASE CF TFE FRECUENCY COMPONENTS. 

C/LL CRANP(NUMPTS*FFEC»FHASE*ITE>fiT8) 

S TC F 
EMC 
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SieRCLTINE WAVE 

IMS SUBROUTINE COMPUTES THE VALUES OF AN Arbitrary 
SELF-GENERATEC COMPLEX WAVEFORM FOR A SINGLE FFECUENCY. 
TFE FREQUENCY IS SPECIFIEC BY THE LSER. 

CCMMCN X(ia24),Y(1024),XW<512l,YW(512>,XS<512l 
CCMMCN YS(512)>XP(512)»YP(51Z)»FREQ(512)»HR(S12) 
CCMMCN WI(512),GR<512),GI<512),ANF<512),TINE<512) 
CCMMCN FHASE<5I2),SS(512)tF,T,N,RN,SPL,FRPV«,TPV.,FCRW 
CCMMCN FCT1»PCT2»FCT3>PCT4»PCT5»PCT6 
I £ = T/M 
CC 30 « = 1,N 
K * M - 1 
J » M * I 

FFEC(J) = M * (l./T) 

X S ( M ) = CCS16.28318 * F * K * T$ ) 

YS<M) = SIM6.2331E * F » K * TSI 
30 CCMTINLE 
RETURN 
EMC 
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SUEFCUTINE PRMULT 

J.^1» ~0BRCLTINE COMPUTES THE X-AXI ' VALUES FCR PLPTTTKr 
CIEM5." Ti C0MPIjTeS TFE PRE-MULTIPLICATI^ CC5FFI- 

CfMMrN «?3!’ vn/ c?2i ’in} 1 12 J*Vte(S12J » XS ( 512 ) 

LCrMCN iS(512)tXP(512) * YP ( 5 1 2 ) • FFECt‘512 ) .UR 1^1? I 

^^°^-^Iij?^2!FH3fpCTi^Fcf5;Fcf6 PL ’ PPPW,TS ‘"’ FaSk 

CCM FLTE TFE TIME INCREMENTS FOP PLOTTING TFE CLTFUT. 
TIME(M) = K * (T/RN) 

CCMFLTE TFE PRE-MULT I FLICAT ION CCE F F IC IENTS . 

S « t S P l J * (K**2)/2. 

C7 F - 6 > F ( S ) 

JMEl = J £12 > * CCS(<3.1415S*(K*=»2) l/P.N) 

VMM) = ( C7R ) * aIN( (3. 1415S*( K**2) )/RNI 



IF CALLEC, TFE 

CALCULATED FOR 

TFE F cFCEN'TAGE OF OEVIATICN FROM A NOMINAL ' VALLE* 



PRS-NLLTIFLICATION COEFFICIENTS ARE 
A UNIFORMLY RANCCM RANGE CETEFM IN EC EY 



1C 



IF ( FF.RW .N E . 1.0) GO TO 
CALL RANCLI IX,IY,YFL) 
? h C “ Y F L 
HIM = XU ( m ) 

YMN) = YVt(N) 

I) = IY 
CCMINLE 
FETIFN 
EK 



i a 



( (pcti)*xw( y )) 

( ( FCT i ) *YW ( M ) ) 



(RNC*PCT2*XM Y) } 
(RNC*eCT2*YMY > ) 
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SLEPOUTINE XSUM 

THE SLERCLT INS MULTIPLIES THE SAMFLEC DATA ev TPS PRE- 
PL IT I PL IC AT I CN COEFFICIENTS ANC SU«S ThE REAL ANC IMAG- 
INARY COMPONENTS. 

CCMMON X ( 1024) * Y( 1024) ,XXI< 512) ,>X( 512), XS (512 1 
CCMMCN >S(512)»XP(512)»YP(512)*FREC(512),XR(512) 
CCMMCN XI (512) ,GR(512),GI<512) , AFP (512) ,TIME(512) 
CCMMCN PHASE ( 512 ) * SS ( 512) ,F ,T ,N ,PN,SPL, PPRX,TR«,PCRX 
CCMMCN FCT1»PCT2,PCT3,PCT4,PCT5,FCT6 
CC 16 P = 1 , N 
X F 1 = XMM!> * XS ( M ) 

X R 2 = XS (M) * (-1. null'll 

>11 * YS(M) * (-1. *Y'MM) > 

>12 = YS (M ) *XMM) 

SLM THE REAL ANO IMAGINARY TERMS CP THE PRE-MU IT I PL I- 
C AT 1 CN PPICP TC CONVOLVING IN TFE CHIRP FILTERS. 

X F ( M I * XR 1 - Y 1 1 
X 1 (M ) = XR 2 + Y I 2 
U CONTINUE 
FETLRN 

enc 



c 

c 

c 

c 



c 

c 

c 

c 

c 



HEFCL7INE TAP* 

This SL6RCLT INE COMPUTES THE TAPPING WEIGHTS FCP THE CO- 
PIE* CCNVCLLTICN FILTERS. 

CCMNCN *(1024) ,Y(1024) ,XW(512) , Y V< < 5 1 2 ) t X $ ( 5 12 ) 

CCNMCN YS(512),XP(512),YP(512)tFPEC(512),WP(512) 

CCNHCN *I(512),GR(512),GI(512) ,AMP(512) ,TIME(512) 
CCNHCN PHASE (512 ) ,SS(512) , F,T, N ,RN ,SPL, PRRW ,TPW,FCP* 
CCMNCN PCTltPCT2,PCT3,PCT4, FCT5 ,FCT6 
FT s N 

CC 11 N = 1,N 
K = N- h 

5 = (($FO* (K**2)/2») * (-1.0) 

CIS = E *P ( S ) 

>(M = (CTR) * CCS((3.14159*(K**2) )/RN) 

Y(N) = (CTR) * SIM(3.1415S*(K**2) )/RN ) 

11 CCMINLE 
N M = N - 1 

CC 14 K = 1,NM1 

2 = N - K 

N = N + K 

> (M = * ( J ) 

Y (M = Y ( J ) 

14 CCNT IM E 

IF CALLEC, THE TAPPING WEIGHTS ARE CALCULATEC FCR A UNI; 
FCP.UY RANCCM 0 ANGE CETERNINEC BY THE PERCENTAGE CF CSVI- 
ATICN FFCH A NCMINAL VALUE. 

IF (TRW.NE.l.O)RETURN 
NPV =2 *N-1 
1> = 111111111 
CC 75 J = 1*NRW 
CALL PANCL(IX»IY f YFL) 

>U') = >(J) - ( (PCT3)*X(J) ) + (*NC1*PCT4**( J) ) 

>(J) = Y(J) - ( ( PCT2 ) * Y ( J ) ) + ( PNC1*PCT4*Y ( J ) ) 

I> = IY 
7 f CCNTINLE 
FETLPN 
EK 
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SIESCLTINE CCNV 

CCNVC LLTKN T INE SHIFTS * MtLTIPLYS, ANC SUMS TC FERFnsp 



rrvUrlu vc » 2 1 3 ! ’ Ik i » l , i } 2 1 * ™ < 512 ) , XS ( 5 12 ) 

C C N NCN YS(512)*XP(512)*YP(512)»FRSC(512) .^R (51 

rr£2r!!! Si 1 ( I J? ] t Gl ( 5 12 ) , AN P ( 512) ,T I V E ( 5 



12 ) 

12 ) 



CC 15 P = 1,N 



HM = C.O 

scya * c.o 

SC! «2 = C.C 

st,»* « = c.o 

CC 20 K = 1, 



I F (N. 


GT 


.1) GO 


TO 


1 2 


J - 


N 


■f 


1 - K 






Gl 




= 


hP (K ) 


* 


X ( J ) 


G 2 
G 2 




s 


WF (K) 


* 


V(J) 




= 


W I ( K ) 




X(J) 


0* 




s 


WI (K) 


* 


V ( J ) 


GC 


TC 


2 










N 


+ 


N - K 






Gl 




= 


W » ( K > 


* 


X < J ) 


G 2 




= 


WR(K) 


* 


V ( J ) 


G 2 




= 


W I ( K ) 


* 


X( J) 


C4 




= 


WI(K) 


* 


Y ( J ) 



'ey LF THE VALUES AFTER SHIFTING ANC NLLTIFLYINC- IN THE 
CCNVC LLTICN FILTERS. 



i SIM 


= Gl 


+ 


SUN 1 


S IN 2 


= G2 


+ 


SUN 2 


SU3 


= G 3 


+ 


SUNS 


SIM 


= C-4 


+ 


SUN 4 



2C CCMINLE 

CCYFLTE THE REAL AND IMAGINARY CCMFCNEM VALLES AT THE 
CLTFLT CF THE CONVOLVER. 

G F < y ) = SLM1 - SUM 4 
Gl <y ) = SLR 2 + SLH3 
If CC N T INL E 
F ETLPN 
ENC 
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JIEFCLTINE POMULT 

JtiKM^SSa'SI^SSIb'SI 

^^- C ^^llIp^2’pCT3% l !l;f^; R F^r UFBRl ’* T,,l '* FCP '' 

ccjs'm'S'ITn 

U Mihir 2,/2 - 

I F ( FCRU .NeIi!o* 



CfU B4fCU(IX,IY,YFL) 

M Li = > F L 

>£ <E! = XP(N,) -hpct5)«xp(m)) 
I> - !Y YP(M) ~ ( <PCT5)*VP(M)) 

CC h 7 I NU E 
PEURN 
E C 



( RfVC2*PCT6 > P ( M | ) 
(RND2*PC76*YP(y) ) 
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S l EFCUT IN E AMPHA 

THIS SUBROUTINE PCST-NLLTIFLIES THE PEAL ANC IMAGINARY 
CCNYCLLTICN FILTER CUTPUTS AND THEN CETERMINES THE M AGNI 
UCE ANC PHASE CF THE FPECUENCY CCMFCNENTS. 

CCMMCN X ( 1024 ) » Y ( 102 A ) , XH ( 5 12 ) , Yfc { 5 12 ) t XS ( 5 1 2 ) 

CCNMCN YS(512),XP(512),YP(512) , PREC ( 512 ),V*P (512) 
CCNMCN WI(512),GR{512)tGI(512>,ANF<512),TIME<512) 
GCNMCN FhASE(512)»SS(512),F,T,N,PN,SPL,PPPV«iTRW,FCRW 
COHGN FCTlt PCT2,PCT3,PCT4, FCT5 ,FC76 
CC SO N = 1 y N 
El = GF (M ) * XW ( M ) 

£2 = GR (M ) * (-1. *YW(M) ) 
c 3 * G 1 ( M ) * (-1. *YW(M) 

ea = g i ( y I *xwm 

SUM THE FCS TNLLT I PLI EC VALUES TC DETAIN THE IMAGINARY ANC 
PEAL VALUES CF THE FPECUENCY COMPONENTS. 



PE 

PI 



a fil - B3 
=82+84 



CETEFMINE THE M AGNITUC6 ANC PHASE CF THE FREQUENCY CC-MPC 
N EM S • 

AFC- = (PE**2) + ( P 1**2 ) 

A N F (M ) = ( (SCRT(APC-) )/PN) 

IF (AES(PE) .LT.O.OCDRE = 0.C 
IF (AES(FI).LT.0.001)RI = 0.0 
FHASE(M) = 100.0 ^ „ 

I FIFE. EC. 0.0. ANC. PI. GE. 0.0) PHASE(M) = 90.0 
IF (RE.EC.C.O. AN D.PI.LT.0.0) PHASE (M ) = -90.0 
I F (P5.G E.Q.O.AND.PI.EQ.O.O) PHASE ( N ) =0.0 
IF (PE.LT.0.0.AN0.P I.EC.C.O) PHASE (M ) = 180. 0 
I F ( PH AS E ( M ) .NE.100.Q ) C-C TC 90 
PHASE (M ) = ATAN2(PI,RE>* 57.2957 
SC CCM1NUE 
FETLPN 
EMC 
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PROGRAM LISTING CF THE PRIME TRANSFORM ALGORITHM 



Tt-IS FROGRAM USES THE "PRIME TRANSFORM” ALGORITHM TC 
CCMRITE THE FOURIER COEFFICIENTS OF SAMPLED PEAL CATA. 

C CM NON >S (32), X( 64 ),V( 64 ),P(32, 32), PT <32, 22), GP (32) 
COMMON NR EM (32) , R c < 22 ) , P I ( 3 2 J , X M A C- ( 32 ) , X ? < 2 2 ) 

C C M MON FR6(32) , PR I (32 ),PHASE(22) ,TIME(32 ),FFEC<32 ) 
CCNMCN M,T,NM,T$ 

INTEGER *4 ITB ( 1 2 ) / 12*0/ 

R E A L*4 RTB ( 28 ) / 28*0 .0/ 

PEAC IN THE CATA WIN'CC’a LENGTH "T M ANC THE NLMEER OF CATA 

samples "N", Inhere "n m is an gcc prime. 

F EAC (5 , 10 )T ,N 

PEAC IN THE R6SICUE VALLES TAKEN ^RCM CRC TABLES FOR 
THE RELATION (R**N) MCC(N). 

NM = N - 1 

P E AC ( 3 , 11) (NREM(K) ,K=1,NM) 

CALL S’aAVE 
CALL FRMLT 
CALL PEFVnAV 
CALL TRANk* 

CALL TPCSE 
CALL XMC-PH 

CLTFLT THE CCMPUTED CATA. 

WFITE (6 ,205 ) 

VnFITE(6 ,240 

WFITE ( 6 , 2 AO ) 4 k 4 

WRITE(6,200)((P( I,J),J=ltNM) ,I=1,NM) 

W F I T 6 ( 6,240) 

Vn F ITE (6 ,215) 

Lp it! ( I lloo! ((PT( I, J) » J = 1 ,NP) , I = 1 , NM ) 
fcfI1E(6 ,24C) 

WfnE(6,260) 
i* s nE < 5 , 22G ) 

V Fill (6 ’,210) ( I,XS( I ),XF (I ) , 1 = 1, N ) 

WFITEI6 , 240 ) 

WfnE<6,235) 

V» F I T E ( 6 ,245) 

VifI7 c (6i?40) 

v!FITE(6’220)(I,NRE)-(n,X(n,Ym,PE(I),^im,FSE(n, 

J F c I Cl ) ,I«l,Nf<) 
l» F J7E 16 » 24C ) 

V* F I T £ (6,255) 

W F ITE ( 6 , 1 10 ) ( I,XMiC-< 1), PHASE ( I) ,I*1,N) 
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I7E (2 ) = 3 
I T E ( 3 ) * 8 
ITE<4) = 5 
I T E 112) = C 
PIE (II = T/8.0 
PTE (2 ) = 2.0 

PICT THE SAMPLED INPUT WAVEFORM. 

C<LL CRAWP(N, TI ME,XS,ITe,RTet 
F T E U ) * 2.0 
N F T = N / 2 
PTE (2 » * 0.2 

PICT THE MAGNITUDE OF THE FRECUENCY CCMPCNENTS. 

CALL C«AWF<NPT,FREC,XMAG, ITE,RTE ) 

«TE(2) = 80.0 

FICT TPE PHASE CF THE FRECUENCY COM FCNENTS . 

K c Fc^^§:^r Ec ’ pHASE ’ iTe ’ RTE ' 

11 PC FN AT { 1615) 

11C FCFMT(/,I5,2F20.3) 

2CC FCFNAT(/,2X,20F4.1) 

2C5 FCFi*AT(40X, ’PERMUTATION MATRIX "F'") 

2 1C FCFNAT(/,2X,I5,3X,2F20.5) 

ill FC FMAT ( 25 X , ' THE TRANSPOSE OF THE FERMUTATICf^ NATPIX') 
PCFFAT 2X,' INDEX’ , 12X, 'NORMAL CP C E FED * , 1 2X , * P ER MO T E C * ) 
2 - C FCFMAT(/,2X,I5,I1C,6F17.3) 

2 2 5 FCFMAT(2X, 'INDEX* , 5X , • RE S I CUE 1 , 5 X , » COS I N 5 V»EIC-HT',5X, 
} 'SINE ‘rjcIul-T 1 ,7X, 'REAL PART', 6 X , * I PAG I NAP V PART', AX, 
1'REAL PART* ,6X, • IMAGINARY PART') 

FCF W AT(//) 

FC FMAT (57X, * PERMUTEC CRCER ' ,4X , • P ERMUTED CRCEF', 

I4X, 'NORMAL ORDER' ,4X, 'N'CRMAL CRCEF') 

FCPMATl 2X, ' INDEX* , 10 X, •MAGNITUDE ' , 12X , ' PHASE* ) 
FCF*AT(24X,' INPUT DATA') 

STOP 
ENC 



24C 

245 

255 
2 5 C 
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SlEFCUTINE S WAVE 



THIS SU3RGIT IN E COMPUTES The SAMPLE VALUES FCR AN ARBI- 
TRARY fcAVEFCRM AS CEFINEC EY ThE USER. IT At$C COMPUTES 
ThE SAMPLE TI^E, THE FUNC AMENT AL FRECUENCY ANC ITS hAR- 
M C N I CS • 

ThE SAMPLE TIMES ANC FReCUENCI ES ARE USEC AS ThE X-AXIS 
IN ThE CLTFUT PLOTS. 

CCMMCN XS(32) ,X(64) , Y ( 64 ) , P (32 , 22 ), PT(32,22 ),GP(22 ) 
CCMMCN NR EM (32) ,RE(32) , RI ( 32 ) , XMAG ( 32 ) , XF ( 3 2 I 
COMCN PR E (32), PR I (32) , PHASE (22) , TIME (32) »F«EC(32 ) 

CCM NON N,T,NM,TS 
TS = T/N 



= 1 , N 



CC 12 M 
K = M - 
TINE(M) 

FFEC(M) 

X S ( M ) = 

I 2 E 212 * 

$ CGS((6. 28318 * 

$- (2.0)*CCS(6.28218 * 
CCNTINUE 
RETURN 
ENC 



= K 
= K 
0.8 
6.0 



TS 

(l./T ) 

COS (6.28218 
K * IS) + 

14. 



* 2. * K * TS ) + ( 2.C ) *S IN ( 6. 

* K * TS) + 0.78529) 



22.0 * K * TS) 
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S l 6 FCUT IN E PRMUT 

THIS SlBRGLTINE CGMPUTES THE PERMUTATION MATFIX «P» THAT 
IS LS6C TC PERMUTE THE SAMPLED REAL CATA. 1HE TRANSPOSE 
CF «P« IS ALSO CGMPITEC. 

CCMMCN XS(32) ,X(64) t Y(6A),P(32,22 )tPT(3 2,22 >tCP(32) 
CCMMCN NREM ( 32 ) » R E ( 3 2 I »RI(32),XMAG(32),XP(22) 

CCMMCN PRE (22 ) ? PRI ( 22 ) , FHAS E ( 2 2 I , T I M £ ( 32 ) , F F E C ( 3 2 ) 

CCMMQN N»T,NM,TS 

CC 20 NN = 1 1 NM 

CC 25 M. = 1 1 NM 

IFINREM (NN).NE.M) GC TC I 

P(NN,M) = 1.0 

GC TC 25 

1 F <NN,M> = C.O 
2 i CCNTINLE 
2C CCNTINLE 



1 P AN S PCS E THE PERMUTED MATRIX 

CC E5 M = 1 * NM 
CC SO K = 1 * NM 
FT (K » M ) = P(M,K) 

SC CCNTINLE 
Sf CCNTINLE 
PETLPN 
ENC 
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S l E RCUT INE PERWAV 

f Ir|i S j|?fLeI NE PES,1,ljTeS THE SANPLEC CATA EXCLLCING THE 

FERf-LTFr pv mF.itt^vt^ * 11 matrix anc is 

re i? MULTIPLYING I< 6Y Thr PEFNUTATICN MTRIX 
UUCP TR 4ft (.\M X NN ( MATRIX. THE RESULTANT NATIXIS 



M-ICP IS AN 
(NN X 1). 



CCNNCN NReM?ii?, ( g|i42>fBli32»f5pl^-2> 3 X^T?zi GP<32) 

EcEScn Kf^if§ RI ‘ 32, - PHASEI32, - TlRE<S2 »'”«‘«> 
£ c i I°cio = l ’ NM 

CC iO J = 1,NM 
K = J + 1 

PC = PC 1 1 J 1, * XS ( K ) 

G 1 = P G + C- 1 
60 CONTINUE 
CP (I) = G 1 
fC CONTINUE 

CC 30 J = 1,NM 
K - J + 1 
XF(1) = XS(1) 

XMK) * GPU) 

20 CCNTINLE 
RETURN 
ENC 
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SIERCLTINE TRANk 

THIS SUBROUTINE COMPUTES THE TAPPING HEIGHTS FCP THE 
COMPLEX CIFCLLAR CORRELATION THAT MAY EE IMPLEMENTED 
USING A TRANSVERSAL FIUTER AND COMPUTES THE RESULTS 
CF THE COR REU AT ION. THE TAPPING HEIGHTS CORRESPOND TC 
THE VALUES OF THE M C" MATRIX. 

C C M MON XS(22),X<64)»Y<64),P<32,32),PT(32,22),GP<22) 
CCMMCN NREV( 32) ,RE(321»RI ( 3 2 ) , XM A G ( 32 I , X F ( 3 2 ) 

COMMON FRE(22) , PR I (32) , PHASE(22) , T IME < 32 ) , F REC ( 32 ) 
COMMON N,T,NN,TS 

CALCULATE THE TRANSVERSAU FILTER CC E F FIC IENTS . 

CC TC K = I » NM 
I = NM 4 K 

X ( K ) - C0S((6. 28318 * NREM ( K ) ) /N ) 



7C 



Y (K ) 

X(I) . . 

Y ( I ) s Y ( K ) 
CONTINUE 



(-1) * SIN ((6. 28319 '* NREM(K))/N) 
> (K ) 



4 M - KK 

* GP(NM+1-KK) 

* GP(NM+1-KK) 



PERMUTED CATA IS CORRELATED IN THE TRANSVERSAL FILTER. 

CC 75 M = 1,NM 
SUM1 * C.C 
SUM2 =• C.Q 
CC 60 NK = 1,NM 

J = NM 4 1 

X F 1 = X(J) 

YF1 = Y(J> 

SUM AFTER SHIFTING ANC MULTIPLYING 

SUM1 = > R 1 4 SUM 1 
SLM2 = YR1 4 SUM 2 
EC CONTINUE 

THE PEAL PART IS OBTAIN EC FROM THE CCSINE FILTER ANC THE 
IMAGINARY PART IS OBTAINED FROM THE SINE FILTER. 

F E (M ) = SLM1 
P I (M) = S L w 2 
75 CCNTINLE 
REURN 
ENC 
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S l E RGLT IN E T-POS E 



7 HIS SUBROUTINE MULTIPLIES THE CIRCULAR CORRELATES RE- 
SULTS ev THE TRANSPOSE OF THE PERMUTATION MATRIX AN C 
THEN A CDS THE FIRST CATA SAMPLE TG EACH CQMFLTcC REAL 
FART. 

C CM MON XS(22),X(64)tY<e4),P<32,22),PT(32t32 ),CP(32J 
CCMMCN NREM(22) ,RE(22) t R I ( 3 2 } , XM A C ( 22 ) , XF(22) 

CCMMCN FRE (32) , PR I (32) , PHASE (22) ,7 XM5( 32 ),FREC( 22 ) 

C C M MON M,T,NM,TS 
CC E5 M = ltNM 
Cl = c.c 
C* 2 * C.C 
CC 100 K = 1 1 NM 
>1 = PT(VtK) * RE ( K ) 

Y 1 = PT ( M t K ) * RI(K) 

Cl = XI ♦ Cl 
G2 = Y1 + C-2 
ICC CCMINLE 

FRE(M) * Cl + XS ( 1 ) 

FFI (M) = G 2 
G5 CCMINLE 
FE1L5N 
EMC 
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SLEFCUTINE XMGPH 

THIS S U6R0 UT I NE COMPUTES THE MAGNITUCE AND PHASE FOR 
EACH FREQUENCY COMPONENT TC GIVE A CNE-SIDEC SPECTRA. 
NOTE THAT THE C-C COMPONENT IS CQMFLTEC SEPARATELY. 

C C M MON XS (32) ,X<64) , Y ( 64 ) , P ( 3 2 , 2 2 ) , PT ( 3 2 , 2 2 ),GP<22 ) 
COMMON NREN<32),RE(32),RK32)tXMAG<22),XP<22) 

COMMON PRE ( 32 ) * PR 1(32) v PHASE ( 32 ) v TIKE ( 32 ) * FPEC ( 32 ) 
COMMON M,T,NM,TS 
CC 115 M = 1,NM 
J = M4l 

IF (AES( PRE(Ml) ).t_T .0.001 ) PRE(MI) = 0.0 
IF ( ABS ( PRI(M) J.LT .0.001 ) PRKM) = C.O 
AFC* = <FRE(M)**2) + <PRKM)**2) 

X M AG ( J ) = ’ (SCRT( ARG) J/N 
FHASE(J) = 100.0 

IF(PRECM).EC.O.O.ANC.PRICM).GE.O.O)PHASECJ) = SO.G 
IF (PRE (M).6Q. 0.0. AND.PRKMJ.LT. C.O PHASE (J ) =-9C.C 
1F(PRE(M). GE. 0.0. AND. PRKM). EO.G.G ) PHASE < J ) = 0.0 
IF (PRE(M).LT. 0.0. ANC. PRKM) .EC. C.OJPHASEU ) = 1S0.C 
IF(FHASE( JJ.NE.1CC.0) C-C TO 115 
PHASE (J ) = ATAN2CFRICK ) ,PRE<M) ) * 57. 2957 
Ilf CONTINUE 

COMPUTE THE "C-C" COM PCNENT MAGNITUCE ANC PHASE. 

PFASE(l) = 0.0 
GC = O.C 
CC 1C 5 . = 1 » N 
C*C = XS ( J ) + GC 
105 CONTINUE 

XMAG(l) = GC/N 

FETURN 

ENC 
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